Python 的欧拉方法

Question

我想用欧拉法在 0 到 2 的区间上逼近 dy/dx = -x +1 的解。我正在使用这段代码

def f(x):
    return -x+1  # insert any function here


x0 = 1  # Initial slope #
dt = 0.1  # time step
T = 2  # ...from...to T
t = np.linspace(0, T, int(T/dt) + 1)  # divide the interval from 0 to 2 by dt
x = np.zeros(len(t))
x[0] = x0 # value at 1 is the initial slope
for i in range(1, len(t)):  # apply euler method
    x[i] = x[i-1] + f(x[i-1])*dt

plt.figure()  # plot the result
plt.plot(t,x,color='blue')
plt.xlabel('t')
plt.ylabel('y(t)')


plt.show()

我可以使用此代码来逼近任何函数在任何区间上的解吗？很难看出这是否真的有效，因为我不知道如何在近似值旁边绘制实际解 (-1/2x^2 + x)。

Answer 1

您可以使用以下方法绘制实际的解决方案：

def F(x):
   return -0.5*x+x

# some code lines

plt.plot(t,x,color='blue')
plt.plot(t,F(t),color='orange')

但请注意，实际解 (-1/2x+x = 1/2x) 与您的斜率 f(x) 并会显示不同的解决方案。

实际解 (-1/2x+x = 1/2x) 的 * 实际斜率 f(x) 是只是 f(x)=1/2

Answer 2

如果您始终对同一角色使用相同的变量名，这可能会有所帮助。根据您的输出，解决方案是 y(t)。因此你的微分方程应该是dy(t)/dt = f(t,y(t))。这将给出斜率函数及其精确解的实现

def f(t,y): return 1-t

def exact_y(t,t0,y0): return y0+0.5*(1-t0)**2-0.5*(1-t)**2

然后将 Euler 循环也作为一个单独的函数来实现，尽可能地保留问题的具体细节

def Eulerint(f,t0,y0,te,dt):
    t = np.arange(t0,te+dt,dt)
    y = np.zeros(len(t))
    y[0] = y0 
    for i in range(1, len(t)):  # apply euler method
        y[i] = y[i-1] + f(t[i-1],y[i-1])*dt
    return t,y

然后将解绘制为

y0,T,dt = 1,2,0.1
t,y = Eulerint(f,0,y0,T,dt)
plt.plot(t,y,color='blue')
plt.plot(t,exact_y(t,0,y0),color='orange')

Python 的欧拉方法

Euler's method for Python

python

matplotlib

differential-equations