# 二阶偏微分方程组 龙格库塔法_备战数模（一）——龙格库塔法解微分方程与编程实现...

## 一、龙格-库塔法解微分方程

Step1:

Step2:将

Step3:

### 3.编程实现

python如下编程：

import numpy as np
import matplotlib.pyplot as plt

class ode_RK_4:
def __init__(self, func, T_span, Y0, Nt):
"""
初始化函数
:param func:需要计算的导函数即f(x,y)
:param T_span: x的范围
:param Y0: Y的初始值
:param Nt: 需要分割的个数
"""
self.func = func
self.T = T_span
self.Y0 = Y0
self.Nt = Nt
self.t_array = np.linspace(T_span[0], T_span[1], Nt)
self.y_array = np.zeros(Nt)
self.y_array[0] = Y0
self.h = (T_span[1] - T_span[0]) / (Nt - 1)

def main(self):
"""
龙格库塔计算函数
:return: 返回x的列表与y的列表
"""
for i in range(0, self.Nt - 1):
K1 = self.func(self.t_array[i], self.y_array[i])
K2 = self.func(self.t_array[i] + self.h / 2, self.y_array[i] + self.h / 2 * K1)
K3 = self.func(self.t_array[i] + self.h / 2, self.y_array[i] + self.h / 2 * K2)
K4 = self.func(self.t_array[i] + self.h, self.y_array[i] + self.h * K3)
self.y_array[i + 1] = self.y_array[i] + self.h / 6 * (K1 + 2 * K2 + 2 * K3 + K4)
return self.t_array, self.y_array

if __name__ == '__main__':
def func(x, y): return np.exp(x) + y

a = ode_RK_4(func, [0, 2], 1, 10)
t, y = a.main()

plt.plot(t, y,lw=2, label="'RK-4's result")
plt.plot(t, np.exp(t) * (t + 1), '-.', label="Real")
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.show()

### 5.高阶方程组具体例子

class Double_P():
"""
这是双复摆的一个类
"""

def __init__(self, m1, m2, l1, l2, g):
A = m1 * l1 ** 2 / 6 + m2 * l1 ** 2 / 2
B = m2 * l2 ** 2 / 6
C = m2 * l1 * l2 / 2
E = l1 * (m1 + 2 * m2) * g / 4
F = l2 * m2 * g / 4
M = np.array([[A / E, C / (2 * (E * F) ** 0.5)], [C / (2 * (E * F) ** 0.5), B / F]])
result = np.linalg.eig(M)
# 下面是一些物理值的计算，在今天的代码中使用不到
lamb1, lamb2 = result[0]
print('lmd1,lmd2')
print(lamb1, lamb2)
omega1, omega2 = 1 / (lamb2) ** 0.5, 1 / (lamb1) ** 0.5
print('omega1,omega2:')
print(omega1, omega2)
print('lambda1,lambda2:')
print((lamb1 - A / E) * 2 * E / C, (lamb2 - A / E) * 2 * E / C)
# 在self中储存一些参量
self.lambda1 = (lamb1 - A / E) * 2 * E / C
self.lambda2 = (lamb2 - A / E) * 2 * E / C
self.omega1 = omega1
self.omega2 = omega2
self.att = [m1, m2, l1, l2, g]

def no_OF(self, t, y):
"""
这是没有外力情况下的纽扣函数
:param t: float
:param y: 当前状态量
:return: ddtheta, ddphi, dtheta, dphi
"""
m1, m2, l1, l2, g = self.att
dtheta, dphi, theta, phi = y
coefficient_matrix = np.array([[m1 * l1 ** 2 / 3 + m2 * l1 ** 2, m2 * l1 * l2 * np.cos(theta - phi) / 2],
[m2 * l1 * l2 * np.cos(theta - phi) / 2, m2 * l2 ** 2 / 3]])
result_matrix = np.array(
[-l1 * (m1 + 2 * m2) * g * np.sin(theta) / 2 - m2 * l1 * l2 * dphi ** 2 * np.sin(theta - phi) / 2
,
-l2 * m2 * g * np.sin(phi) / 2 + m2 * l1 * l2 * dtheta ** 2 * np.sin(theta - phi) / 2])
ddtheta, ddphi = line_s.solve(coefficient_matrix, result_matrix)
return ddtheta, ddphi, dtheta, dphi

def sol_ode(self, t_max, y0, step, func, args=None, return_Q=False):
"""
this is func to cal double_bai's theta and phi with time going
:param t_max: the max_t
:param y0: ini_attribute:dtheta0, dphi0, theta0, phi0
:param step: time_step in ode
:param attribute: m1,m2,l1,l2,g
:return: a array like array_t np.arange(0, t_max, step)
"""
t_span = [0, t_max]
t_eval = np.arange(0, t_max, step)
result = sci_int.solve_ivp(func, t_span=t_span, y0=y0, max_step=0.01, t_eval=t_eval,
method='DOP853', args=args)
return result.y[:, :], t_eval

01-11

06-01
03-26
02-11
01-06
01-16 1058
10-28 5542
10-17 1万+
08-22 1万+