integral function returns 0 value
3 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
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 comentarios
Dyuman Joshi
el 11 de Oct. de 2023
@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.
Respuestas (1)
Walter Roberson
el 12 de Oct. de 2023
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)
min(Kvec_sym), double(ans)
max(Kvec), double(ans)
max(Kvec_sym), double(ans)
0 comentarios
Ver también
Categorías
Más información sobre Calculus en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!