计算科学 - Matlab 中的简单集成已关闭 - 吾爱随笔录

Matlab 中的简单集成已关闭

计算科学 matlab 一体化

2021-11-28 18:38:23

这是一个令人难以置信的平凡集成，但由于某种原因，我无法发现错误。我只是在球坐标中执行一个简单的 3d 积分，它应该返回 1.0 的值：

Nr = N;
Nt = round(pi*N);
Np = round(2*pi*N);

rs = linspace(R0, Rf, Nr);
ts = linspace(0, pi, Nt);
ps = linspace(0, 2*pi, Np);

dr = rs(2)-rs(1);
dt = ts(2)-ts(1);
dp = ps(2)-ps(1);

%fprintf('Total number of nodes: %i\n', Nr*Nt*Np);

C = 1/((4/3)*pi);

fint = 0.0;
for ir = 2:Nr
  r = rs(ir);
  r2dr = r*r*dr;
  for it = 1:Nt-1
    t = ts(it);
    sintdt = sin(t)*dt;
    for ip = 1:Np-1
      p = ps(ip);
      fint = fint + C*r2dr*sintdt*dp;
    end 
  end 
end

这与

\int_{0}^{2 π} \int_{0}^{π} \int_{0}^{1} C r^{2} \sin (θ) d r d θ d ϕ

$\int_0^{2\pi} \int_0^\pi \int_0^1 C r^2 \sin(\theta) dr\, d\theta\, d\phi$

我在哪里

C = \frac{1}{(4 / 3) π}

$C = \frac{1}{(4/3)\pi}$ 因为这是球体体积的倒数。

但是，我发现我返回的整数值有一个一致的错误，大约是(1/N)/(1.0 - numeric_integral) ~= 0.66. 所以我在某个地方犯了一个错误，但我不知道在哪里或为什么。

1个回答

您正在使用一阶精确积分技术（矩形规则）并且您的误差与 $1/N\propto\Delta r,\Delta \phi,\Delta \theta$ . 这正是您应该期望的收敛行为类型。

如果您使用更准确的积分方案，您将看到更快的收敛。您可以轻松地为这个积分实现梯形规则，我希望看到二阶收敛： $\text{Error}\propto(1/N)^2$ .

编辑以添加更多详细信息：

您使用的矩形规则与以下相同（检查您的答案，它们应该与舍入误差相匹配）：

%create a 3d grid of values
[PHI,THETA,R] = ndgrid(ps,ts,rs);

%function to be integrated
F = C*R.^2.*sin(THETA);
fint = sum(sum(sum(F(1:end-1,1:end-1,2:end))))*dr*dt*dp;

而梯形规则可以实现为

%create a 3d grid of values
[PHI,THETA,R] = ndgrid(ps,ts,rs);

%function to be integrated
F = C*R.^2.*sin(THETA);

% INTEGRATE:
% r direction
I = sum( F(:,:,1:end-1)+F(:,:,2:end), 3 );
% theta direction
I = sum( I(:,1:end-1)+I(:,2:end), 2 );
% phi direction
I = sum( I(1:end-1)+I(2:end) );
fint = I*dr*dp*dt/8;

计算不同的误差 $N$ 并在对数图上绘制揭示了它们的准确性顺序。正如所料，梯形规则是 $\mathcal{O}(1/N^2)$ 而矩形规则是 $\mathcal{O}(1/N)$ .

对数误差图

其它你可能感兴趣的问题

上一篇在时间积分器中准确计算当前时间下一篇Nesterov 的第一种分段线性凸优化问题方法中的 beta