用python语言实现三阶和四阶龙格库塔法求解方程
时间: 2023-11-25 14:07:53 浏览: 68
好的,下面是用 Python 实现三阶和四阶龙格库塔法求解方程的代码:
```python
import numpy as np
def f(t, y):
"""
定义方程f(t, y)
"""
return t * np.sqrt(y)
def rk3(t0, y0, h, T):
"""
三阶龙格库塔法求解方程
"""
N = int(T / h) # 步数
t = np.linspace(t0, T, N+1)
y = np.zeros(N+1)
y[0] = y0
for i in range(N):
k1 = h * f(t[i], y[i])
k2 = h * f(t[i] + h / 2, y[i] + k1 / 2)
k3 = h * f(t[i] + 3 * h / 4, y[i] + 3 * k2 / 4)
y[i+1] = y[i] + (2 * k1 + 3 * k2 + 4 * k3) / 9
return t, y
def rk4(t0, y0, h, T):
"""
四阶龙格库塔法求解方程
"""
N = int(T / h) # 步数
t = np.linspace(t0, T, N+1)
y = np.zeros(N+1)
y[0] = y0
for i in range(N):
k1 = h * f(t[i], y[i])
k2 = h * f(t[i] + h / 2, y[i] + k1 / 2)
k3 = h * f(t[i] + h / 2, y[i] + k2 / 2)
k4 = h * f(t[i] + h, y[i] + k3)
y[i+1] = y[i] + (k1 + 2 * k2 + 2 * k3 + k4) / 6
return t, y
```
其中,`f(t, y)` 函数定义了方程的右侧,`rk3(t0, y0, h, T)` 函数用三阶龙格库塔法求解方程,`rk4(t0, y0, h, T)` 函数用四阶龙格库塔法求解方程。这两个函数都接受四个参数:
- `t0`:初始时刻;
- `y0`:初始值;
- `h`:步长;
- `T`:结束时刻。
调用这两个函数可以得到数值解。例如,要求解方程 $y' = t\sqrt{y}$,初始条件为 $y(0)=1$,在 $[0,1]$ 区间内,步长为 $0.1$,可以这样调用:
```python
t3, y3 = rk3(0, 1, 0.1, 1)
t4, y4 = rk4(0, 1, 0.1, 1)
```
最后,我们可以使用 Matplotlib 将数值解绘制成图像:
```python
import matplotlib.pyplot as plt
plt.plot(t3, y3, label='RK3')
plt.plot(t4, y4, label='RK4')
plt.legend()
plt.xlabel('t')
plt.ylabel('y')
plt.show()
```
相关推荐
![zip](https://img-home.csdnimg.cn/images/20210720083736.png)
![zip](https://img-home.csdnimg.cn/images/20210720083736.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)