integral function returns 0 value
3 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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 commentaires
Dyuman Joshi
le 11 Oct 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.
Réponses (1)
Walter Roberson
le 12 Oct 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 commentaires
Voir également
Catégories
En savoir plus sur Calculus dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!