simple problem with double integral, why give me an error?

6 views (last 30 days)
hello, I should solve a double integral in which the integrand function contains parameters that I must subsequently obtain through a minimization.
f=
the code I wrote is the following:
clc
clear all
x0=[1 1 1 1 1 1]
lb=[0.01 0.01 0.01 0.01 0.01 0.01]
ub=[5 5 5 5 5 5 ]
ftomin=@(x) provaInt(x(1),x(2),x(3),x(4),x(5),x(6))
options = optimoptions('simulannealbnd','FunctionTolerance', 1e-6,'MaxIterations',500,'MaxFunEvals',500,'PlotFcns',{@saplotbestx,@saplotbestf,@saplotx,@saplotf,@saplottemperature,@saplotstopping});
[x,fval,exitflag,output]=simulannealbnd(ftomin,x0,lb,ub,options)
function f=provaInt(St,at,bt,am,bm,Sm)
syms t m
fi=((1./(St.*(t-3)*sqrt(2*pi))).*exp((-1/2).*((log(t-3)-at-bt*3)./St).^2)).*((1/(sqrt(2*pi)*Sm)).*exp((-1/2).*((m-am-bm*4)./(Sm)).^2)));
f=int(int(fi,t,1970,1980),m,5.75,8));
%where all parameters to be found are positive and real (not imaginary numbers)
%f=-(2535301200456458802993406410752*pi*(erf((2^(1/2)*(1/Sm^2)^(1/2)*(am + 4*bm - 8))/2) - erf((2^(1/2)*(1/Sm^2)^(1/2)*(am + 4*bm - 23/4))/2))*(erfi((2^(1/2)*(-1/St^2)^(1/2)*(at + 3*bt - log(1967)))/2) - erfi((2^(1/2)*(-1/St^2)^(1/2)*(at + 3*bt - log(1977)))/2)))/(31859534503965572279823959492121*Sm*St*(1/Sm^2)^(1/2)*(-1/St^2)^(1/2))
end
All the parameters to be obtained with the optimization are positive and real (no imaginary numbers).
If I supply the integral with the "int" function (as in the example) matlab gives me the following error:
Conversion to logical from sym is not possible.
Error in acceptancesa (line 31)
if delE < 0
Error in sanewpoint (line 30)
if options.AcceptanceFcn(optimvalues,newx,newfval)
Error in saengine (line 29)
solverData = sanewpoint(solverData,problem,options);
Error in simulanneal (line 50)
solverData = saengine(solverData,problem,options);
Error in simulannealbnd (line 163)
[x, fval, exitflag, output] = simulanneal(FUN, x0, [], [], [], [], lb, ub, options);
Error in prova_funzione_integrale (line 11)
[x,fval,exitflag,output]=simulannealbnd(ftomin,x0,lb,ub,options)
Caused by:
Failure in AnnealingFcn evaluation.
If, on the other hand, I supply the integral already resolved ("f" commented in green) the optimization takes place without problems.
In reality I have to solve a quadruple integral in which the extremes of integration are defined and are numbers. However, to better understand how to do this type of calculations I tried with a double integral with a part of the function that I have to solve.
What am I doing wrong?
Thank you very much for helping.

Accepted Answer

Steven Lord
Steven Lord on 28 Sep 2020
The function you pass into simulannealbnd must return a double array not a sym array. Convert the value of the integral into double using double before returning it from your function.
  3 Comments
Steven Lord
Steven Lord on 28 Sep 2020
You could try solving the integral once symbolically and use matlabFunction to turn the result into a function that computes the answer numerically, not symbolically. However the symbolic solution may be more complicated or MATLAB may not be able to generate the answer without you telling it that it can assume certain things about the variables. For example, if I asked you to draw the curve y = a*x^2 you can't do that without me telling you the sign of a so you know whether it opens up or down (or is the straight line y = 0.)
Emanuele Biondini
Emanuele Biondini on 29 Sep 2020
Thank you so much for your help. I tried to fix the program to try to speed it up. However, the computation time, although improved, is still excessively slow. In the example the "for" loop that is used to make a summation, goes from 1 to 237, in the "true" program that I have to set up this number is much larger, for example 100'000. Is there a way to speed up the computation time?
This is the fixed code
Thank you
clc
clear all
load('Dataset_Aquila.mat')
mi=Ctemp(:,10);
ti=Ctemp(:,14);
xi=Ctemp(:,8);
yi=Ctemp(:,7);
L=length(mi);
x0=[0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 ]
lb=[0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01]
ub=[9 9 9 9 9 9 9 9]
ftomin=@(x) provaInt(x(1),x(2),x(3),x(4),x(5),x(6),x(7),x(8),mi,ti,xi,yi)
options = optimoptions('simulannealbnd','FunctionTolerance', 1e-6,'MaxIterations',500,'MaxFunEvals',500,'PlotFcns',{@saplotbestx,@saplotbestf,@saplotx,@saplotf,@saplottemperature,@saplotstopping});
[x,fval,exitflag,output]=simulannealbnd(ftomin,x0,lb,ub,options)
function fD=provaInt(St,at,bt,am,bm,Sm,Sa,ba,mi,ti,xi,yi)
lo=log(10);
sqr2pi=sqrt(2*pi);
Sum=0;
for i=1:237
i
fi=@(x,y,t,m)((lo.*exp((-lo.*(am+(bm-1)*mi(i)+Sm.^2.*lo/2)))).*((((t-ti(i))./(St.*(t-ti(i))*sqr2pi)).*exp((-1/2).*((log(t-ti(i))-at-bt*mi(i))./St).^2))).*((1/(sqr2pi*Sm)).*exp((-1/2).*((m-am-bm*mi(i))./(Sm)).^2))).*((1./(2.*pi.*(Sa.^2).*(10.^(ba*mi(i)))))).*exp(-((x-xi(i)).^2+(y-yi(i)).^2)./(2.*(Sa.^2).*10.^(ba.*mi(i))));
%iterm=integral(@(m)integral3(@(x,y,t)fi(x,y,t,m),min(xi),max(xi),min(yi),max(yi),2000,2020),5.75,8,'ArrayValued',true)
iterm=integralN(fi,min(xi),max(xi),min(yi),max(yi),2000,2020,5.75,8)
tf = isreal(iterm)
if tf==0
disp('error')
pause
end
Sum=Sum+iterm;
end
fD=vpa(-Sum)
end

Sign in to comment.

More Answers (2)

Walter Roberson
Walter Roberson on 29 Sep 2020
I did not do a detailed analysis, but at the moment it looks to me as if iterm can only be complex valued if
log(t-ti(i))
is complex valued, which would be the case ti(i) > t. If I am correct then you could test
if ti(i) > 2000
and fail the situation immediately without doing the integral() step.
Also,
Sum=0;
That is double precision
iterm=integralN(fi,min(xi),max(xi),min(yi),max(yi),2000,2020,5.75,8)
integralN calculates in double precision
Sum=Sum+iterm;
so that will be double precision
fD=vpa(-Sum)
and there you invoke the symbolic toolbox to calculate the rational approximation to the double precision number and return that (unless the double precision number happens to be "approximately" a square root of an interesting number, in which case it will return a square-root form.)
ftomin=@(x) provaInt(x(1),x(2),x(3),x(4),x(5),x(6),x(7),x(8),mi,ti,xi,yi)
so your objective function is going to do its calculation in floating point, but at the last minute it converts it to symbolic.
[x,fval,exitflag,output]=simulannealbnd(ftomin,x0,lb,ub,options)
to a routine that operates in double precision and will promptly convert the rational approximation into double precision number -- one that might be slightly different than the negative of the double precision sum calculated by the function.
The trip into symbolic and back out again wastes time and reduces precision. Just have your function return -Sum instead of vpa(-Sum)
  2 Comments
Emanuele Biondini
Emanuele Biondini on 29 Sep 2020
Thanks Walter for your kind reply. No, on the contrary, t-t (i) must be positive. Complex numbers must not come out.
iterm=integralN(fi,min(xi),max(xi),min(yi),max(yi),2000,2020,5.75,8)
tf = isreal(iterm) %possible solution: 1=iterm is real, 0=iterm is imaginary and the program stops
if tf==0
pause
end
The part where I evaluate if the number is real, I wrote it just for verification.
Maybe I could try to solve the integral first leaving the result in parametric form, in this way it doesn't have to be recalculated every time.
Walter Roberson
Walter Roberson on 29 Sep 2020
No, on the contrary, t-t (i) must be positive
It would fail to be positive if t < ti(i) which is what I expressed as "which would be the case ti(i) > t" -- that this, in context, the failure would happen in that case. So your "on the contrary" does not seem correct, as we are expressing the same thing ?

Sign in to comment.


Emanuele Biondini
Emanuele Biondini on 29 Sep 2020
Sorry, I misinterpreted the sentence. Yes, we are saying the same thing. I deleted vpa (-sum) but the computation time is still very very long. I also apologize if I am slow to understand, but I had never done this type of calculations with matlab..

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!