清空筛选条件
清空筛选条件

integral function returns 0 value

3 次查看（过去 30 天）
park junho
park junho 2023-10-11
回答： Walter Roberson 2023-10-12
Ran in:
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 个评论
显示 1更早的评论
Dyuman Joshi
Dyuman Joshi 2023-10-11
@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
park junho 2023-10-12
Ran in:
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 个）

Walter Roberson
Walter Roberson 2023-10-12
Ran in:
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

类别

在 Help Center 和 File Exchange 中查找有关 Calculus 的更多信息

标签

产品


版本

R2020a

Community Treasure Hunt

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

Start Hunting!

Translated by