如何在 MATLAB 中对无限极限进行数值积分?
How to numerically integrate with infinite limit in MATLAB?
我想对一个无穷大的积分进行数值积分。有谁知道我应该怎么做?
int(x* exp (v*x + (1-exp(v*x))/v),x, o , inf)
无效。
请注意,我将具有 v
的值。
%n=10;
kappa=.5;
delta0=.5;
Vmax=500;
Vdep=2.2;
l=2.2;
kbT=4.1;
%xb=.4;
fb=10;
k=1;
V0=5;
e1=(fb*l/kbT)*(kappa/delta0);
e2=Vmax/V0;
e3=Vdep/V0;
w=zeros(1,25);
for v=1:25
w(:,v)=integral(@(x) x.*exp(v*x+((1-exp(v*x))/v)),0,inf);
end
e12=e2*exp(-e1*(1:25).*w.^2)-e3;
plot(e12);
ylim([0 25]);
hold on;
plot(0:25,0:25);
xlim([0 25]);
%hold off;
该图与文章中的真实数据不符!(特别针对e12曲线)
我需要计算两条曲线的交点(根据论文为 ~13.8)然后在第二部分我必须在 e12 中添加一个包含自变量的项:
v=13.8;
w= integral(@(x) x.*exp(v*x+((1-exp(v*x))/v)),0,inf)
e4 = zeros (1,180);
fl = 1:180;
e4(:,fl)= (fl*l/kbT)*(kappa/n);
e12=e2*exp(-e1*v*w^2-e4)-e3
但同样的问题是 运行 这段代码我将以 e12 的负值结束,它应该在大的 fl (fl>160)
中接近零
要显示此代码与预期曲线有何不同,您可以将这些数据绘制在同一张图上:
fl = [0, 1, 4, 9, 15, 20, 25, 40, 60, 80, 100, 120, 140, 160, 180];
e12 = [66, 60, 50, 40, 30, 25.5, 20, 15.5, 10.5, 8.3, 6.6, 5, 2.25, 1.1, 0.5];
这显然与代码生成的曲线不符。
通过对距离为 dx
的离散点处的函数求和来执行数值积分。您选择的 dx
越小,您获得的近似值就越好。例如,从 x=0
到 x=10
的积分是通过以下方式完成的:
x = 0:dx:10;
I = sum(x.* exp (v*x + (1-exp(v*x))/v))*dx;
显然,您不能为 x=inf
这样做。但我相信你的功能会迅速衰退。因此,您可以假设 x* exp (v*x + (1-exp(v*x))/v) = 0
足够大 x
。否则积分发散。所以你所要做的就是设置 x
的限制。如果您不确定限制应该是多少,您可以执行一个带有停止条件的循环:
I = 0;
prevI = -1;
x = 0;
while abs(I-prevI)>err
prevI = I;
I = I + x.* exp (v*x + (1-exp(v*x))/v)*dx;
x = x + dx;
end
现在,您只需设置所需的 dx
和 err
你必须阅读这篇文章:Mathwork Link
也许您使用的函数有误。另请注意,MATLAB 语法区分大小写..
假设问题是关于这个完整代码的:
syms x;
v = 1; % For example
int(x*exp(v*x + (1-exp(v*x))/v),x, 0, Inf)
问题是它 returns 本身(即 int
没有找到解析解),可以设置 'IgnoreAnalyticConstraints'
option to true
(more details) 来得到一个解决方案:
syms x;
v = 1; % For example
int(x*exp(v*x + (1-exp(v*x))/v),x, 0, Inf, 'IgnoreAnalyticConstraints', true)
returns -ei(-1)*exp(1)
,其中 ei
is the exponential integral function (see also expint
for numerical calculations). For negative values of v
the solution will also be in terms of eulergamma
, the Euler-Mascheroni constant。当然,如果 v
是 0
.
,积分是未定义的
使用 Mathematica 10.0.2 的 Integrate
生成符号 v
.
的完整解决方案
Integrate[x Exp[v x - (Exp[v x] - 1)/v], {x, 0, Infinity}]
returns
ConditionalExpression[(E^(1/v) (EulerGamma + Gamma[0, 1/v] + Log[1/v]))/v, Re[v] < 0]
正在应用Assumptions
:
Integrate[x Exp[v x - (Exp[v x] - 1)/v], {x, 0, Infinity}, Assumptions -> v > 0]
Integrate[x Exp[v x - (Exp[v x] - 1)/v], {x, 0, Infinity}, Assumptions -> v < 0]
returns
(E^(1/v) Gamma[0, 1/v])/v
和
(E^(1/v) (2 EulerGamma - 2 ExpIntegralEi[-(1/v)] + Log[1/v^2]))/(2 v)
其中 Gamma
is the upper incomplete gamma function。这些看起来与 Matlab 的结果相匹配。
要在 Matlab 中对这些进行数值计算:
% For v > 0
v_inv = 1./v;
exp(v_inv).*expint(v_inv).*v_inv
或
% For v < 0
v_inv = 1./v;
exp(v_inv).*(2*double(eulergamma)+2*(expint(v_inv)+pi*1i)+log(v_inv.^2)).*v_inv/2
我想对一个无穷大的积分进行数值积分。有谁知道我应该怎么做?
int(x* exp (v*x + (1-exp(v*x))/v),x, o , inf)
无效。
请注意,我将具有 v
的值。
%n=10;
kappa=.5;
delta0=.5;
Vmax=500;
Vdep=2.2;
l=2.2;
kbT=4.1;
%xb=.4;
fb=10;
k=1;
V0=5;
e1=(fb*l/kbT)*(kappa/delta0);
e2=Vmax/V0;
e3=Vdep/V0;
w=zeros(1,25);
for v=1:25
w(:,v)=integral(@(x) x.*exp(v*x+((1-exp(v*x))/v)),0,inf);
end
e12=e2*exp(-e1*(1:25).*w.^2)-e3;
plot(e12);
ylim([0 25]);
hold on;
plot(0:25,0:25);
xlim([0 25]);
%hold off;
该图与文章中的真实数据不符!(特别针对e12曲线) 我需要计算两条曲线的交点(根据论文为 ~13.8)然后在第二部分我必须在 e12 中添加一个包含自变量的项:
v=13.8;
w= integral(@(x) x.*exp(v*x+((1-exp(v*x))/v)),0,inf)
e4 = zeros (1,180);
fl = 1:180;
e4(:,fl)= (fl*l/kbT)*(kappa/n);
e12=e2*exp(-e1*v*w^2-e4)-e3
但同样的问题是 运行 这段代码我将以 e12 的负值结束,它应该在大的 fl (fl>160)
中接近零要显示此代码与预期曲线有何不同,您可以将这些数据绘制在同一张图上:
fl = [0, 1, 4, 9, 15, 20, 25, 40, 60, 80, 100, 120, 140, 160, 180];
e12 = [66, 60, 50, 40, 30, 25.5, 20, 15.5, 10.5, 8.3, 6.6, 5, 2.25, 1.1, 0.5];
这显然与代码生成的曲线不符。
通过对距离为 dx
的离散点处的函数求和来执行数值积分。您选择的 dx
越小,您获得的近似值就越好。例如,从 x=0
到 x=10
的积分是通过以下方式完成的:
x = 0:dx:10;
I = sum(x.* exp (v*x + (1-exp(v*x))/v))*dx;
显然,您不能为 x=inf
这样做。但我相信你的功能会迅速衰退。因此,您可以假设 x* exp (v*x + (1-exp(v*x))/v) = 0
足够大 x
。否则积分发散。所以你所要做的就是设置 x
的限制。如果您不确定限制应该是多少,您可以执行一个带有停止条件的循环:
I = 0;
prevI = -1;
x = 0;
while abs(I-prevI)>err
prevI = I;
I = I + x.* exp (v*x + (1-exp(v*x))/v)*dx;
x = x + dx;
end
现在,您只需设置所需的 dx
和 err
你必须阅读这篇文章:Mathwork Link
也许您使用的函数有误。另请注意,MATLAB 语法区分大小写..
假设问题是关于这个完整代码的:
syms x;
v = 1; % For example
int(x*exp(v*x + (1-exp(v*x))/v),x, 0, Inf)
问题是它 returns 本身(即 int
没有找到解析解),可以设置 'IgnoreAnalyticConstraints'
option to true
(more details) 来得到一个解决方案:
syms x;
v = 1; % For example
int(x*exp(v*x + (1-exp(v*x))/v),x, 0, Inf, 'IgnoreAnalyticConstraints', true)
returns -ei(-1)*exp(1)
,其中 ei
is the exponential integral function (see also expint
for numerical calculations). For negative values of v
the solution will also be in terms of eulergamma
, the Euler-Mascheroni constant。当然,如果 v
是 0
.
使用 Mathematica 10.0.2 的 Integrate
生成符号 v
.
Integrate[x Exp[v x - (Exp[v x] - 1)/v], {x, 0, Infinity}]
returns
ConditionalExpression[(E^(1/v) (EulerGamma + Gamma[0, 1/v] + Log[1/v]))/v, Re[v] < 0]
正在应用Assumptions
:
Integrate[x Exp[v x - (Exp[v x] - 1)/v], {x, 0, Infinity}, Assumptions -> v > 0]
Integrate[x Exp[v x - (Exp[v x] - 1)/v], {x, 0, Infinity}, Assumptions -> v < 0]
returns
(E^(1/v) Gamma[0, 1/v])/v
和
(E^(1/v) (2 EulerGamma - 2 ExpIntegralEi[-(1/v)] + Log[1/v^2]))/(2 v)
其中 Gamma
is the upper incomplete gamma function。这些看起来与 Matlab 的结果相匹配。
要在 Matlab 中对这些进行数值计算:
% For v > 0
v_inv = 1./v;
exp(v_inv).*expint(v_inv).*v_inv
或
% For v < 0
v_inv = 1./v;
exp(v_inv).*(2*double(eulergamma)+2*(expint(v_inv)+pi*1i)+log(v_inv.^2)).*v_inv/2