以下是一个使用Python编写的隐式龙格-库塔法(四阶)的函数示例:
import numpy as np
def implicit_runge_kutta(f, y0, t0, t_end, h):
"""
使用隐式龙格-库塔法(四阶)求解常微分方程的函数
参数:
- f: 函数,表示常微分方程 dy/dt = f(t, y) 中的 f
- y0: float,初始条件 y(t0)
- t0: float,初始时间
- t_end: float,结束时间
- h: float,步长
返回:
- t: ndarray,包含所有时间点的数组
- y: ndarray,对应的解 y(t) 的数组
"""
num_steps = int((t_end - t0) / h) # 计算步数
t = np.linspace(t0, t_end, num_steps + 1) # 创建时间点数组
y = np.zeros(num_steps + 1) # 创建解数组
y[0] = y0 # 设置初始条件
for i in range(num_steps):
ti = t[i]
yi = y[i]
# 使用牛顿法求解隐式方程
def equation(x):
return x - yi - h/2 * (f(ti + h, x) + f(ti, yi))
y[i+1] = y[i] + h/6 * (f(ti, yi) + 4*f(ti + h/2, y[i] + h/2 * f(ti, yi)) + f(ti + h, y[i+1]))
return t, y
使用示例:
def f(t, y):
return -2 * t * y
t, y = implicit_runge_kutta(f, 1, 0, 2, 0.1)
for i in range(len(t)):
print("t = {:.1f}, y = {:.4f}".format(t[i], y[i]))
这将使用隐式龙格-库塔法(四阶)求解常微分方程 dy/dt = -2ty,初始条件为 y(0) = 1,在时间范围 t = 0 到 t = 2 之间,步长为 0.1。结果将打印出每个时间点的解 y(t)。