在python 3中使用RK4方法求解二维线性流的耦合微分方程

计算科学 计算物理学 软件 微分方程 线性系统
2021-12-27 22:28:01

我想研究二维线性流的动力学,其动力学方程是-(x1˙x2˙)=(1142)(x1x2). 现在,我尝试在 python 中使用 RK4 作为初始条件来求解并绘制这个耦合微分方程的 y 与 x 的关系(y0=2,x0=1). 我的代码如下,但图形不正确[在图形原点应该是鞍点,(0.251)轴应该是稳定的歧管和(11)轴应该是不稳定的歧管]-

import numpy as np
from math import sqrt

import matplotlib.pyplot as plt



# Equations:
   
def V(u,t):
    x1, x2, v1, v2 = u
    return np.array([ v1, v2, (x1+x2), -(4*x1-2*x2)])

def rk4(f, u0, t0, tf , n):
    t = np.linspace(t0, tf, n+1)
    u = np.array((n+1)*[u0])
    h = t[1]-t[0]
    for i in range(n):
        k1 = h * f(u[i], t[i])    
        k2 = h * f(u[i] + 0.5 * k1, t[i] + 0.5*h)
        k3 = h * f(u[i] + 0.5 * k2, t[i] + 0.5*h)
        k4 = h * f(u[i] + k3, t[i] + h)
        u[i+1] = u[i] + (k1 + 2*(k2 + k3) + k4) / 6
    return u, t


u, t  = rk4(V, np.array([2.0, 0., -1, 1.]) , 0. , 1. , 1000)
x1, x2, v1, v2 = u.T
plt.plot(x1,x2)
plt.show()

谁能帮我写这段代码。

1个回答

您正在实施系统

(x¨1x¨2)=(1142)(x1x2)
在其一阶版本中vk=x˙k. 这种不同的系统导致不同的解决方案也就不足为奇了。

对于您声明的系统,您必须使用

def V(u,t):
    x1, x2 = u
    return np.array([  (x1+x2), (4*x1-2*x2)])

状态空间维度的差异应该很明显。