integral function returns 0 value

Hello.
I would like to fit the transient decay data with the multi exponential function convoluted with the gaussian irf function.
before using lsqcurvefit, I checked my integration function work well but it return 0 value.
if I change the xdata_fit date to xdata_fit = linspace(100,1500,500), it return correct results.
What is wrong in my code and how to handle it?
xdata_fit = linspace(100,1500,100);
parameter = [0.0,0.3,9,0.18,5,300];
y0=parameter(1);
t0=parameter(2);
A=parameter(3);
t1=parameter(4);
B=parameter(5);
t2=parameter(6);
f3 = @(x) (A.*(1-exp(-x/t1))+B.*exp(-x/t2)).*exp((-(xdata_fit-x-t0).^2)*4*log(2)/(0.15^2));
Kvec = integral( @(x) f3(x),0,inf,'Arrayvalued',true)+y0; % Original
figure()
plot(xdata_fit,Kvec)

3 个评论

Can you imagine how small the expression
exp(-(xdata_fit-x-t0).^2)
becomes when xdata_fit starts at 100 ?
@Torsten I believe the function values can't be represented in double precision, so they are effectively 0.
And so, the integrals are (effectively) 0 as well.
Thanks for comments.
exp(-(xdata_fit-x-t0).^2) gonna be very small in x=100 to 1500.
However, as I said in question, if I set the x data to xdata_fit = linspace(100,1500,500), integration is working well!
xdata_fit = linspace(100,1500,500);
parameter = [0.0,0.3,9,0.18,5,300];
y0=parameter(1);
t0=parameter(2);
A=parameter(3);
t1=parameter(4);
B=parameter(5);
t2=parameter(6);
f3 = @(x) (A.*(1-exp(-x/t1))+B.*exp(-x/t2)).*exp((-(xdata_fit-x-t0).^2)*4*log(2)/(0.15^2));
Kvec = integral( @(x) f3(x),0,inf,'Arrayvalued',true)+y0; % Original
figure()
plot(xdata_fit,Kvec)

请先登录,再进行评论。

回答(1 个)

The numeric integration is too low of a precision and is giving a result that is very wrong.
format long g
Q = @(v) sym(v);
xdata_fit = linspace(100,1500,500);
parameter = [0.0,0.3,9,0.18,5,300];
y0=parameter(1);
t0=parameter(2);
A=parameter(3);
t1=parameter(4);
B=parameter(5);
t2=parameter(6);
f3 = @(x) (A.*(1-exp(-x/t1))+B.*exp(-x/t2)).*exp((-(xdata_fit-x-t0).^2)*4*log(2)/(0.15^2));
Kvec = integral( @(x) f3(x),0,inf,'Arrayvalued',true)+y0; % Original
%switch to symbolic
parameter = Q(parameter);
y0=parameter(1);
t0=parameter(2);
A=parameter(3);
t1=parameter(4);
B=parameter(5);
t2=parameter(6);
syms x
f3sym(x) = (A.*(1-exp(-x/t1))+B.*exp(-x/t2)).*exp((-(xdata_fit-x-t0).^2)*4*log(Q(2))/(Q(0.15)^2));
Kvec_sym = vpaintegral(f3, x, 0, inf) + y0;
Kvec_symd = double(Kvec_sym);
figure()
plot(xdata_fit,Kvec); title('original numeric');
figure()
plot(xdata_fit,Kvec_symd); title('symbolic');
min(Kvec), double(ans)
ans =
1.44241509999502
ans =
1.44241509999502
min(Kvec_sym), double(ans)
ans = 
6.0713773021437683471964604310212e-97159106
ans =
0
max(Kvec), double(ans)
ans =
2.009645779592
ans =
2.009645779592
max(Kvec_sym), double(ans)
ans = 
0.0000000000000000000000037657495170743369911765134795177
ans =
3.76574951707434e-24

类别

帮助中心File Exchange 中查找有关 Creating and Concatenating Matrices 的更多信息

产品

版本

R2020a

标签

Community Treasure Hunt

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

Start Hunting!

Translated by