Python06 向前Euler法、向后Euler法、梯形方法、改进的Euler方法以及四阶Runge

1 解如下常微分方程

2 分别使用向前 Euler 法、向后 Euler 法、梯形方法、改进的 Euler 方法以及 四阶 Runge_Kutta 方法 结果如下图所示

由结果可以发现 向前 Euler 法、向后 Euler 法和准确值的误差比较大 梯形方法、改进的 Euler 方法以及四阶 Runge_Kutta 方法和准确值更为接近 并 且在该微分方程中的效果相差不大。
3 因为上图中梯形方法、改进的 Euler 方法以及四阶 Runge_Kutta 方法的重叠 单独的结果如下图

from class6 import ode
ode ode() 
ode.shows()

import math
import numpy as np
import matplotlib.pyplot as plt
class ode(object):
 包括: 1. 显示Euler; 2. 隐式Euler; 3. 梯形方法; 4. 改进Euler方法; 5. Runge-Kutta方法(四阶) 
 def __init__(self):
 # 初始值
 self.y0 1
 # 区间[a,b]
 self.a 0
 self.b 1
 # 步长
 self.h 0.1
 # 定义常微分方程
 def ODE(self, x, y):
 f -y x 1
 return f
 # 计算准确值
 def acc(self, x):
 # 方程真解为 y exp(-x) x
 acc math.exp(-x) x
 return acc
 # 1. 显示Euler法
 def Euler_forward(self):
 yn self.y0
 xn self.a
 accuracy []
 Euler_forward []
 x []
 while xn 1:
 Euler_forward.append(yn)
 accuracy.append(self.acc(xn)) 
 x.append(xn)
 f self.ODE(xn, yn)
 # 显示Euler方法: y(n 1) yn hf(x(n 1),y(n 1)) , x(n 1) x0 (n 1)*h x(n) h
 yn yn self.h * f
 xn xn self.h
 return x, Euler_forward, accuracy
 # 2. 隐式Euler
 def Euler_back(self):
 yn self.y0
 xn self.a
 Euler_back []
 x []
 while xn 1:
 xn1 xn self.h
 # 隐式Euler方法: y(n 1) yn hf(x(n 1),y(n 1)) , x(i) x0 i*h xn h
 # 由公式计算可得: y(n 1) (h*x(n 1) yn h)/(1 h)
 yn (self.h*xn1 yn self.h) / (1 self.h)
 Euler_back.append(yn)
 x.append(xn1)
 xn xn1
 return x, Euler_back
 # 3. 梯形方法
 def trapezoid(self):
 yn self.y0
 xn self.a
 trapezoid []
 x []
 while xn 1:
 xn1 xn self.h
 f self.ODE(xn, yn)
 ynn yn self.h * f
 ynnn (self.h*xn1 yn self.h) / (1 self.h)
 # 梯形公式: 向前Euler法和向后Euler方法做算术平均
 yn (ynn ynnn) / 2
 trapezoid.append(yn)
 x.append(xn1)
 xn xn1
 return x, trapezoid
 # 4. 改进Euler方法
 def impr_Euler(self):
 yn self.y0
 xn self.a
 impr []
 x []
 while xn 1:
 xn1 xn self.h
 f self.ODE(xn, yn)
 ynn yn self.h * f 
 f1 self.ODE(xn1, ynn)
 # 改进的Euler公式: 用显示公式作预估 用梯形公式作校正
 yn yn self.h*(f f1) / 2
 impr.append(yn)
 x.append(xn1)
 xn xn1
 return x, impr
 # 5. Runge-Kutta方法(四阶)
 def R_K(self):
 yn self.y0
 xn self.a
 R_K []
 x []
 while xn 1:
 R_K.append(yn)
 x.append(xn)
 k1 self.ODE(xn, yn)
 k2 self.ODE(xn, yn) self.h/2 - self.h*k1/2
 k3 self.ODE(xn, yn) self.h/2 - self.h*k2/2
 k4 self.ODE(xn, yn) self.h - self.h*k3
 # 四阶R_K公式: y(n 1) yn h*(k1 2*k2 2*k3 k4)/6
 yn yn self.h * (k1 2*k2 2*k3 k4) / 6
 xn xn self.h
 return x, R_K