【Python算法】数值分析—四阶荣格库塔方法

4399资深专家

已于 2022-05-04 22:24:03 修改

阅读量3.7k

点赞数 8

分类专栏：数学物理编程文章标签： python pycharm 算法

于 2022-04-07 17:11:26 首次发布

本文链接：https://blog.csdn.net/qq_50920297/article/details/124020783

版权

数学物理编程专栏收录该内容

11 篇文章 15 订阅

订阅专栏

一、前言

四阶荣格库塔方法是一种常微分方程的数值解法，给定区间，划分次数，初值条件后，可以求解出任一点的原函数值。

给定常微分方程初值问题：

$\frac{dy}{dx} = f(x, y), \quad a\le x\le b$

$y(a)=a,\quad h=\frac{b-a}{N}$

利用四阶荣格库塔方法：

$K_{1} =f(x_{n},y_{n})$

$K_{2} =f(x_{n}+\frac{h}{2} ,y_{n}+\frac{hK_{1} }{2} )$

$K_{3} =f(x_{n}+\frac{h}{2} ,y_{n}+\frac{hK_{2} }{2} )$

$K_{4} =f(x_{n}+h,y_{n}+hK_{3} )$

$y_{n+1}=y_{n}+\frac{h}{6} (K_{1} +2K_{2} +2K_{3} +K_{4} )$

可逐次求出微分方程初值问题的数值解 $y_{n}, \quad n=1,2,\dots ,N$

二、代码实现

2.1计算步长及区间数组

def get_list_x_and_h(a,  # 区间起始点
                     b,  # 区间结束点
                     n   # n是划分的次数
                     ):
    """ 获取步长及计算区间 """
    list_x = []
    h = abs(a - b) / n
    for i in list(range(0, n + 1)):
        list_x.append(a + i * h)
    return list_x, h  # h是步长

2.2计算K值

def get_k(fx,  # 方程表达式
          x,
          y,
          precision_value=5  # 计算结果保留的有效数字
          ):
    """ 计算k值 """
    return fx.evalf(subs={'x': x, 'y': y}, n=precision_value)

2.3主函数，循环迭代

def method_runge_kutta(fx,  # 方程表达式
                       list_fx,  # 存储fx值的list
                       a, b, n):
    """ 主循环进行计算 """
    list_x, h = get_list_x_and_h(a, b, n)
    i = 0
    while i < n:
        k1 = get_k(fx, list_x[i], list_fx[i])
        k2 = get_k(fx, list_x[i] + h / 2, list_fx[i] + h / 2 * k1)
        k3 = get_k(fx, list_x[i] + h / 2, list_fx[i] + h / 2 * k2)
        k4 = get_k(fx, list_x[i] + h, list_fx[i] + h * k3)

        value_fx = list_fx[i] + h / 6 * (k1 + 2 * k2 + 2 * k3 + k4)
        list_fx.append(value_fx)
        i += 1

2.4初始化配置

# 设置方程表达式
x, y = sp.symbols('x, y')
fx = sp.sin(x) + y
# 设置迭代初始值
init_fx = -1
list_fx = [init_fx]
# 调用四阶荣格库塔方法
method_runge_kutta(fx=fx,
                   list_fx=list_fx,
                   a=0,
                   b=1,
                   n=5)
# 简单打印迭代结果
print("方程：", fx)
print("迭代结果：", list_fx)

测试结果：

完整代码

import sympy as sp


def get_list_x_and_h(a,  # 区间起始点
                     b,  # 区间结束点
                     n   # n是划分的次数
                     ):
    """ 获取步长及计算区间 """
    list_x = []
    h = abs(a - b) / n
    for i in list(range(0, n + 1)):
        list_x.append(a + i * h)
    return list_x, h  # h是步长


def get_k(fx,  # 方程表达式
          x,
          y,
          precision_value=5):  # 计算结果保留的有效数字
    """ 计算k值 """
    return fx.evalf(subs={'x': x, 'y': y}, n=precision_value)


def method_runge_kutta(fx,  # 方程表达式
                       list_fx,  # 存储fx值的list
                       a, b, n):
    """ 主循环进行计算 """
    list_x, h = get_list_x_and_h(a, b, n)
    i = 0
    while i < n:
        k1 = get_k(fx, list_x[i], list_fx[i])
        k2 = get_k(fx, list_x[i] + h / 2, list_fx[i] + h / 2 * k1)
        k3 = get_k(fx, list_x[i] + h / 2, list_fx[i] + h / 2 * k2)
        k4 = get_k(fx, list_x[i] + h, list_fx[i] + h * k3)

        value_fx = list_fx[i] + h / 6 * (k1 + 2 * k2 + 2 * k3 + k4)
        list_fx.append(value_fx)
        i += 1


# 设置方程表达式
x, y = sp.symbols('x, y')
fx = sp.sin(x) + y
# 设置迭代初始值
init_fx = -1
list_fx = [init_fx]
# 调用四阶荣格库塔方法
method_runge_kutta(fx=fx,
                   list_fx=list_fx,
                   a=0,
                   b=1,
                   n=5)
# 简单打印迭代结果
print("方程：", fx)
print("迭代结果：", list_fx)

如有任何疑问，请在评论区留言，我会努力回复。

源码github地址：https://github.com/method_runge_kutta

gitee地址：https://gitee.com/darlingxyz/method_runge_kutta/tree/master