为什么指数积分器方法的结果不佳？

Question

//u' + Au = g(t,u) can be solved by exponential integrators also
//Following snippet is for exp INtegrators
A = -full(Strang(11)) 
A[end,1]=1;A[1,end]=1;
g(t,u) = 2-u
u0 = zeros(11);u0[6]=1
nsteps = 1000
tmax = 10.0
h = tmax/nsteps
u = u0
t = 0
for k in 1:nsteps
    u = expm(-h*A)*u + h*((expm(-h*A)-1)\(-h*A))*g(t,u)
    t = k*h
end
//this is for euler's method
for k in 1:nsteps
    u += h*(A*u + h*g(t,u))
    t = k*h
end

为什么他们给出的结果很差？

这个方法爆炸的很厉害，应该会收敛到[1.99]*11之类的吧？

执行Exp Integrator时有没有错误？

Answer 1

exponential Euler-Rosenbrock step should be 为 u'=Lu+g(t,u)

u = expm(h*L)*u + ((expm(h*L)-1)/L)*g(t,u)

注意倒矩阵是L（在Matlab中是A/B == A*inv(B)和A\B == inv(A)*B），公式中没有裸露的h。在你的情况下，L = -A.

拆分指数部分的另一种方法如下。设u(t) = exp(-t*A)*v(t) 转化为微分方程

exp(-t*A)*v'(t) = g(t,exp(-t*A)*v((t))

现在应用 v

的正向欧拉步骤

v(t+h) = v(t) + h * exp(t*A)*g(t,exp(-t*A)*v((t))

翻译回 u 中的术语给出

u(t+h) = exp(-(t+h)*A)*v(t+h) 
       = exp(-h*A) * ( u(t) + h * g(t,u(t)) )

这同样会产生一阶方法。

Answer 2

测试题为奇异矩阵。更好的测试是设置：

using SpecialMatrices
A = -full(Strang(11))
g(t,u) = 2-u
u = zeros(11);u[6]=1
nsteps = 10000
tmax = 1.0
h = tmax/nsteps
t = 0

使用这个，修复欧拉中的 h 得到（注意有一个额外的 h，我的错误：

u = zeros(11);u[6]=1
for k in 1:nsteps
    u += h*(A*u + g(t,u))
    t = k*h
end
@show u

u = [0.93573,1.19361,1.26091,1.29627,1.34313,1.37767,1.34313,1.29627,1.26091,1.19361,0.93573]

但要找出问题所在，请先查看数字。 A=0 会发生什么？好吧，我们知道 phi(z) = (e^z - 1)/z。根据 L'Hopital 的规则，phi(z) -> 1 为 z->0。因此，为了使我们的实现具有相同的行为，我们必须具有相同的结果。但是让我们看看会发生什么：

expm(zeros(5,5))

5×5 Array{Float64,2}:
 1.0  0.0  0.0  0.0  0.0
 0.0  1.0  0.0  0.0  0.0
 0.0  0.0  1.0  0.0  0.0
 0.0  0.0  0.0  1.0  0.0
 0.0  0.0  0.0  0.0  1.0

注意这给出了单位矩阵。所以想一想极限：如果底部趋于零……这怎么可能是恒定的？我们必须让顶部变为零...所以顶部将变为 I.

这就是清晰的时刻：作者的意思是 1 在您所在的领域。因此对于矩阵输入，1=I。当你意识到这一点时，你修复了代码：

# Norsett-Euler
u = zeros(11);u[6]=1
for k in 1:nsteps
  u = expm(h*A)*u + ((expm(h*A)-I)/A)*g(t,u)
  t = k*h
end
@show u

u = [0.935722,1.1936,1.26091,1.29627,1.34313,1.37768,1.34313,1.29627,1.26091,1.1936,0.935722]

故事的寓意：对于数学编程，您还必须调试数学。

编辑

一步一步获得更高效的形式。首先，尝试强制使用另一个 varphi 项：

# Norsett-Euler
u = zeros(11);u[6]=1
for k in 1:nsteps
  u = (I + A*(expm(h*A)-I)/A)*u + ((expm(h*A)-I)/A)*g(t,u)
  t = k*h
end
@show u

现在集合：

# Norsett-Euler
u = zeros(11);u[6]=1
for k in 1:nsteps
  u = u + ((expm(h*A)-I)/A)*(A*u + g(t,u))
  t = k*h
end
@show u

这是您尝试编写的方法的有效形式。然后你可以缓存运算符，因为 A 是常量：

# Norsett-Euler
u = zeros(11);u[6]=1
phi1 = ((expm(h*A)-I)/A)
for k in 1:nsteps
  u = u + phi1*(A*u + g(t,u))
  t = k*h
end
@show u

为什么指数积分器方法的结果不佳？

Why is exponential integrator method giving poor results?

matlab

octave

julia

differential-equations

runge-kutta

为什么他们给出的结果很差？

这个方法爆炸的很厉害，应该会收敛到[1.99]*11之类的吧？

编辑