solve nonlinear equations with integration where the upper limit is variable

Hi All,

I need to solve a group of nonlinear equations with integration where the upper limit is variable. The format is shown as below:

f(x,t) is a nonlinear equation without an explicit integraton form. f(x,t) is:

where a1-a6 is known. a6 is more than 800 so it makes f(x,t) have no explicit integraton form.

x(1), x(2), ...x(N) and Z(1), Z(2), ...Z(N) are known double values. The upper limit of the integration, d, is a non-negative variable vector and also the one I need to calculate from the equations. I need to find the smallest ||d|| which satisfies these equations, and || || means any norm of vection d. I plan to use fsolve to solve these equations but I encounter problems in writing the right form of these equations into matlab function handle. My idea is to build the single integration function by int and use for loop to get a N*N matrix M where the element is one single integration. Then I can sum the rows of this N*N matrix and add Z into it. By converting the final function to function handle (matlabFunction) and I can use fsolve. My current diffculty lies in writing the corrrect form of the single integration function. Can anyone give me a hint or tell me is there better idea to solve these equations?

Here is what i arrived

L=10;
interval=1;
x=[-L:interval:L-interval]; %make a x vector between 0 and 10 with 1 as interval
N=length(x);
Z=x.^2./10+0.5;
X=zeros(N);% xi-xj matrix
for i=1:N
    for j=1:N
        X(i,j)=x(i)-x(j);
    end
end
fun=@(dw,xx) integral(@(t)-(a(1)+a(2).*exp(-t./a(3))).*(1-(xx).^2./(a(4)+a(5)*t).^2).^a(6),0, dw, 'ArrayValued',false);
d = sym('d',[1 N]);
D=repmat(d,N,1);% make int upper limit matrix
M = sym('M',[N N]);
for i=1:N
    for j=1:N
        M(i,j)=fun(D(i,j),X(i,j));
    end
end
fun1=sum(M,2)+Z;
fun2=matlabFunction(fun1,'Vars',{[d(1),d(2),d(3),d(4),d(5),d(6),d(7),d(8),d(9),d(10),d(11),d(12),d(13),d(14),d(15),d(16),d(17),d(18),d(19),d(20)]});
d_0=zeros(N,1);
y=fsolve(fun2,d_0)

[John D'Errico]

Shudder.

It would help if you tell people what the function f really is. For example, many people think that a Gaussian hs no known integrable form, and they would be almost correct, in a sense. However, it is entirely tractable, because the integral of a Gaussian exists in the form of the special function erf.

So if you are willing to explain more about the form of f, you might be surprised. For example, if f was a spline, then again, the answer is trivial.

Next, you say you are using int to do the integrations. Of course, that must fail if f has no integrable form.

If you want help, then you need to make it possible to help you. Otherwise, break the problem don into parts. Solve each part. Then put it all together. Don't worry about crap like arrayfun if that is too much for you. Just make it work, and only then worry about making it work more elegantly, more efficiently. This is a common mistake. Elegant crap is meaningless and worthless if you personally don't know how to make it work. So break it all apart. Use loops where you need to use them.

[john zhang]

Hi John, many thanks for your prompt reply and valuable suggestions. I have given the form of my function f and shown my current codes. Hope it helps people to know my problem.

[John D'Errico]

My first comment is you are using symbolic variables in a place where there is no need.

Just because you don't know the value of the vector elements in d, does not mean it must be symbolic. That forces you to jump hurdles, making a symbolic vector, then converting everything back into a function with matlabFunction. You are doing a lot of work that you don't need to do. And worse, because these are syms, it goes MUCH more slowly.

You are using integral and fsolve, not int or solve/vpasolve anyway. For example, in a function like this:

F = @(x) x.^2;

you do not know the value of x that minimizes the function. So what? You could still minimize it as a function of x. You can evaluate the function. You can compute an integral. But I never need to make anything symbolic. Or, suppose I want to find the upper limit of integration such that the definite integral is 17?

Fint = @(u) integral(F,0,u);

I can evaluate that function. (In fact, I can even make it work in a vectorized form if I need to.)

Fint(pi)
ans =
           10.335498066432

The upper limit, where the integral is 17 is then:

fzero(@(u) Fint(u) - 17,10)
ans =
          3.70842976926619

I never needed to make anything symbolic. Of course the answer is correct, since in this case I know the answer.

nthroot(3*17,3)
ans =
          3.70842976926619

But the point is, I never needed to make anything symbolic. Don't confuse not knowing something with needing it to be symbolic.

[john zhang]

Thank you John. According to your hint and suggestion, I descard the use of symbolic variable and try to break the problem into small parts. Now I plan to calculate the first element by using fsolve. Here is my code:

clc
clear

a=[0.0025
    0.0715
    2.4271
   38.8184
   -0.2403
  828.9618];

L=10;
interval=1;
x=[-L:interval:L-interval]; %make a x vector between 0 and 10 with 1 as interval
N=length(x);
Z=x.^2./10+0.5;
X=zeros(N);% xi-xj matrix
for i=1:N
    for j=1:N
        X(i,j)=x(i)-x(j);
    end
end
Fints1=@(d)0;
for i=1:N
    F=@(t)-(a(1)+a(2).*exp(-t./a(3))).*(1-(X(1:i)).^2./(a(4)+a(5)*t).^2).^a(6);
    Fint=@(d)integral(F,0,d(i));
    Fints1=@(d)Fints1(d)+Fint(d(i));
end
d_0=zeros(N,1);
ds=fsolve(Fints1,d_0)

But I encounter this error:

"Error in fsolve (line 230)

fuser = feval(funfcn{3},x,varargin{:});

Error in vtif_test1 (line 33)

dw=fsolve(Fints1,d_0)

Caused by:

Failure in initial objective function evaluation. FSOLVE

cannot continue."

Can you possibly give me some hints? Thank you in advance!

[John D'Errico]

It is difficult to follow what you are doing, and understand what you want to do. For example, I see very early in your code:

x=[-L:interval:L-interval]; %make a x vector between 0 and 10 with 1 as interval

This in fact creates a vector that runs from -10 to +9, with a stride of 1. As such, completely inconsistent with your comment. I'll suppose it creates a vector, I can ignore why you did what you did.

If you have a recent release, thus R2016b or later, the matrix X can be simply created as just:

X = x - x.';

Again, that is irrelevant. But beyond that, what you did, and what you seem to be asking in your question look to diverge rather wildly.

As I try to understand your code and your question, looking back and reading what you have written in your question now leads me to a strange spot, where you are raising numbers to the power of a6, where a6 is 828.9618. And depending on the values of the variables inside, this will result in either underflows or overflows, thus effectively either 0, inf or -inf. Even though you were originally using symbolic variables for part of this, at this level you are still using double precision arithmetic because you are calling integral. And it will with very high probability fail to produce any meaningful result from this computation.

So at this point in your journey, I'll wish you good luck and godspeed, as you will need both in your quest.

[john zhang]

Hi John, thanks for your attention and time. The x=[-L:interval:L-interval] is not that important but "X = x - x.'" does help my coding.

Now I can solve the equations by a most primative method, i.e. adding all the integral and write them in function handle and use fsolve. Though laboriou and time-consuming, I did get one reasonable solution.

Thank you again for your time and attention!