integral function returns 0 value
3 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
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 Kommentare
Dyuman Joshi
am 11 Okt. 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.
Antworten (1)
Walter Roberson
am 12 Okt. 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 Kommentare
Siehe auch
Kategorien
Mehr zu Calculus finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!