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)

[Torsten]

Can you imagine how small the expression

exp(-(xdata_fit-x-t0).^2)

becomes when xdata_fit starts at 100 ?

[Dyuman Joshi]

@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.

[park junho]

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)

Answer 1

[Walter Roberson]

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