Python的Euler方法

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

如果您对同一个角色始终使用相同的变量名，这可能会有所帮助。根据您的输出，解决方案是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')

Answer 2

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

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