链接中的问题: 可以进行解析积分,答案是 4,但是我对使用 Matlab 对它进行数值积分很感兴趣,因为它在形式上类似于我无法进行解析积分的问题。数值积分的困难是因为两个内积分中的函数是
x
、y
的函数,且无法分解出z
和z
。
这绝不是优雅的。希望有人能比我更好地利用matlab函数。我尝试过暴力方式来练习数值积分。我试图通过利用 z=0 处的内部积分也乘以 z 的事实来避免出现极点。我得到 3.9993。必须有人通过使用比梯形规则更好的东西来获得更好的解决方案
function []=sofn
clear all
global x y z xx yy zz dx dy
dx=0.05;
x=0:dx:1;
dy=0.002;
dz=0.002;
y=0:dy:1;
z=0:dz:2;
xx=length(x);
yy=length(y);
zz=length(z);
s1=0;
for i=1:zz-1
s1=s1+0.5*dz*(z(i+1)*exp(inte1(z(i+1)))+z(i)*exp(inte1(z(i))));
end
s1
end
function s2=inte1(localz)
global y yy dy
if localz==0
s2=0;
else
s2=0;
for j=1:yy-1
s2=s2+0.5*dy*(inte2(y(j),localz)+inte2(y(j+1),localz));
end
end
end
function s3=inte2(localy,localz)
global x xx dx
s3=0;
for k=1:xx-1
s3=s3+0.5*dx*(2/(localy+localz));
end
end
嗯,这很奇怪,因为在发帖者之前类似的问题上,我声称这是不可能完成的,现在在看了 Guddu 的答案后,我意识到它并不那么复杂。我之前写过,数值积分会产生一个数字,而不是一个函数,这是正确的——但不是重点:我们可以定义一个函数来计算每个给定参数的积分,这样就可以有效地确实有数值积分结果的函数。
无论如何,就这样吧:
function q = outer
f = @(z) (z .* exp(inner(z)));
q = quad(f, eps, 2);
end
function qs = inner(zs)
% compute \int_0^1 1 / (y + z) dy for given z
qs = nan(size(zs));
for i = 1 : numel(zs)
z = zs(i);
f = @(y) (1 ./ (y + z));
qs(i) = quad(f, 0 , 1);
end
end
我应用了我自己在评论中建议的简化,消除了 x。函数
inner
计算 y 上的内积分值作为 z 的函数。然后函数outer 计算z 上的外积分。我通过让积分从 eps
而不是 0 运行来避免 z = 0 处的极点。结果是
4.00000013663955
inner
必须使用 for
循环来实现,因为赋予 quad
的函数需要能够同时返回多个参数值的值。