欧拉方法的数值稳定性答案

【问题标题】：Numerical stability of Euler's Method欧拉方法的数值稳定性
【发布时间】：2021-04-28 16:09:18
【问题描述】：

我正在尝试估计以下方程组何时在数值上稳定：

dx/dt = -x + y
dy/dt = -100*y

在这个不等式中求解 h

|1+h*λ| <=1. where λ=-1,-100

（方程组的特征值）。我得到 1/50>= h>= 0。欧拉解应该是数值稳定的。然而，当我正在绘制时间步长 h=1/50 的数值和解析解。他们似乎不同意（而对于 h=1e-3，结果看起来不错）。难道我做错了什么？这是我正在使用的代码：

size=20
import matplotlib.pyplot as plt

# Define generl odel

w=1

def Euler(dt,y,x):
        
    dx = dt*(-x+y)
    dy     = -100*y*dt  
    y      = y+dy
    x  = x +dx
    return x,y


# Solve Diff Ep.

tstop =10.0
dt=1/50

y=np.zeros(int(tstop/dt))
x=np.zeros(int(tstop/dt))
t=np.zeros(int(tstop/dt))

t[0]     = 0
x[0]     = 0.1
y[0]     = 0.1



for i in range(0,int(tstop/dt)-1):
        
    dx     = dt*(-x[i]+y[i])
    dy     = -100*y[i]*dt  
    y[i+1]   = y[i]+dy
    x[i+1]   = x[i] +dx
    t[i+1]   = t[i]+dt
   

ax=plt.subplot(1,2,1)
plt.plot(t1,s[:,0],'g-', linewidth=3.0,label='X, Analytical solution')
plt.plot(t,x,"r--",label='X, Euler Method')
plt.plot(t1,s[:,1],'black',linewidth=3.0,label='Y, Analytical solution')
plt.plot(t,y,"b--",label='Y, Euler Method')
plt.text(0.2,0.1, r'$ {\Delta}t \ = \ $'+str(round(dt,4))+'s',horizontalalignment='center',
     verticalalignment='center', transform = ax.transAxes, fontsize=size)

ax.tick_params(axis='x',labelsize=17)
ax.tick_params(axis='y',labelsize=17)
plt.ylabel('X', fontsize=size)   
plt.xlabel('t [sec]', fontsize=size) 
plt.ylim([0,2])
plt.legend(frameon=False,loc=1,fontsize=18)
plt.xscale('log')


plt.ylim([0,0.12])

ax1.tick_params(axis='x',labelsize=17)
ax1.tick_params(axis='y',labelsize=17)
plt.rcParams["figure.figsize"] = [8,8]

【问题讨论】：

我投票结束这个问题，因为它不是编程问题，可能更适合数学 SE。

标签： python numerical-methods ode

【解决方案1】：

使用λ=-100 和h=1/50，您可以获得y 分量的欧拉方法的传播因子值1+h*λ=-1。请注意，解决方案在振荡时保持有界。这就是 A 稳定性的定义。要使收敛为零，您需要h<0.02。要获得非振荡行为，您需要h<=0.01。要查看稍微弯曲的行为，您需要 h<=0.005。

【讨论】：