Python 给定初值，二阶微分方程使用 Runge-kutta 方法

最新推荐文章于 2023-05-02 17:31:51 发布

IcePeak8

最新推荐文章于 2023-05-02 17:31:51 发布

阅读量3k

点赞数

分类专栏：数值分析文章标签： python

本文链接：https://blog.csdn.net/BF2233/article/details/106023164

版权

数值分析专栏收录该内容

1 篇文章

订阅专栏

原函数为
$\begin{aligned} &x''(t)+4x'(t)+5x(t)=0, \\&x(0)=3,x'(0)=-5\\ &\frac{dx}{dt}=y,\\ &\frac{dy}{dt}=-5x-4y \end{aligned}$

import math
import matplotlib.pyplot as plt
import numpy as np
def rks4(F,a,b,Za,M):
    h=(b-a)/M
    T=np.zeros(M+1)
    Z=np.zeros((M+1,len(Za)))
    T=np.linspace(a,b,num=M+1)
    Z[0,:]=Za
    for i in range(M):
        k1=F(T[i],Z[i,:])       
        k2=F(T[i]+h/2,Z[i,:]+h*k1/2)
        k3=F(T[i]+h/2,Z[i,:]+h*k2/2)
        k4=F(T[i]+h,Z[i,:]+h*k3)       
        Z[i+1,:]=Z[i,:]+h*(k1+2*k2+2*k3+k4)/6      
    return T, Z
def F(t,Z):
    x=Z[0]
    y=Z[1]
    return np.array([y, -5*x-4*y])
T,Z=rks4(F,0.,5,[3,-5],50)
actual=np.zeros_like(T)
for i in range(50):
    t=T[i]
    actual[i]=3*math.exp(-2*t)*math.cos(t)+math.exp(-2*t)*math.sin(t)
plt.plot(T,Z[:,0],'r')
plt.plot(T,actual,'b')