为什么指数积分器方法的结果很差？

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

测试问题是一个奇异矩阵。更好的测试是设置：

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

Answer 2

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)) )

这应该同样产生一阶方法。