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

3237ee922d7b8b51b6aa9765ea58ad9b.png

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

在比赛中,我们经常会遇到不可解的微分方程,此时对微分方程求数值解就很有必要。因此,介绍微分方程的一般解法——龙格库塔法。

1.基本思想

考察以下格式,其中

为步长。

希望能够确定系数

使在
前提下使得:

Step1:

做泰勒展开:

Step2:将

带入第一步得:

Step3:

点处的泰勒展开比较:

带入

式得到:

有三个方程,四个个未知量,因此可以根据情况确定。

以上介绍的是龙格库塔的基本思想,下面介绍常用的四阶龙格库塔法

2.四阶龙格-库塔法

类似的,使用不同点值进行组合得到

最常用的是四阶经典龙格库塔法

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()

得到的结果如下:

可以看出结果是很精确地。

c7cfa7c8fa47e5647052c947c5852a0e.png

4.一阶方程组与高阶方程组数值解

方程与高阶方程只需要考虑向量形式就好:

记为

如果太过于复杂,我建议使用库函数:scipy.integrate.solve_ivp 函数官方文档

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

解得画图结果:

640ba646f55d28e96361a3efb4be0371.png
表情包
插入表情
评论将由博主筛选后显示,对所有人可见 | 还能输入1000个字符
相关推荐
©️2020 CSDN 皮肤主题: 1024 设计师:白松林 返回首页