二阶初值问题常微分方程的求解

import numpy as np
from matplotlib import pyplot as plt
plt.rcParams['font.sans-serif']=['SimHei']
plt.rcParams['axes.unicode_minus'] = False


def myfun1(x, y1, y2):
    return y2


def myfun2(x, y1, y2):
    return 2*y2 - 4*y1 + np.cos(x)


def ode1(myfun1, myfun2, x, y1, y2, h, b):
    y1_r = []
    y2_r = []
    x_r = []
    y1_r.append(y1)
    y2_r.append(y2)
    x_r.append(x)
    for i in range(int(b / h)):
        x_old = x
        y1_old = y1
        y2_old = y2
        y1 = y1 + h*myfun1(x_old, y1_old, y2_old)
        y2 = y2 + h*myfun2(x_old, y1_old, y2_old)
        x = x + h
        print('x={:3.2f}\ty1={:3.8f}'.format(x, y1))
        print('x={:3.2f}\ty2={:3.8f}'.format(x, y2))
        y1_r.append(y1)
        y2_r.append(y2)
        x_r.append(x)
    return x_r, y1_r, y2_r


def ode2(myfun1, myfun2, x, y1, y2, h, b):
    y1_r = []
    y2_r = []
    x_r = []
    y1_r.append(y1)
    y2_r.append(y2)
    x_r.append(x)
    for i in range(int(b / h)):
        k11 = h * myfun1(x, y1, y2)
        k12 = h * myfun2(x, y1, y2)
        k21 = h * myfun1(x + h / 2, y1 + 0.5 * k11, y2 + 0.5 * k12)
        k22 = h * myfun2(x + h / 2, y1 + 0.5 * k11, y2 + 0.5 * k12)
        k31 = h * myfun1(x + h / 2, y1 + 0.5 * k21, y2 + 0.5 * k22)
        k32 = h * myfun2(x + h / 2, y1 + 0.5 * k21, y2 + 0.5 * k22)
        k41 = h * myfun1(x + h, y1 + 0.5 * k31, y2 + k32)
        k42 = h * myfun2(x + h, y1 + 0.5 * k31, y2 + k32)
        y1 = y1 + (k11 + 2*k21 + 2*k31 + k41)/6
        y2 = y2 + (k12 + 2*k22 + 2*k32 + k42)/6
        x = x + h
        print('x={:3.2f}\ty1={:3.8f}\tk11={:5.7f}\tk21={:5.7f}\tk31={:5.7f}\tk41={:5.7f}'
              .format(x, y1, k11, k21, k31, k41))
        print('x={:3.2f}\ty2={:3.8f}\tk12={:5.7f}\tk22={:5.7f}\tk32={:5.7f}\tk42={:5.7f}'
              .format(x, y2, k12, k22, k32, k42))
        y1_r.append(y1)
        y2_r.append(y2)
        x_r.append(x)
    return x_r, y1_r, y2_r


if __name__ == "__main__":
    print('欧拉法数值计算结果：')
    x_r, y1_r, y2_r = ode1(myfun1, myfun2, 0, 0, -0.2, 0.1, 1)
    print('标准的四级四阶龙格-库塔法数值计算结果：')
    x_r1, y1_r1, y2_r1 = ode2(myfun1, myfun2, 0, 0, -0.2, 0.1, 1)
    plt.plot(x_r, y1_r, '2:', markersize=1.5)
    plt.plot(x_r, y2_r, '3:', markersize=1.5)
    plt.plot(x_r1, y1_r1, '4:', markersize=1.5)
    plt.plot(x_r1, y2_r1, '8:', markersize=1.5)
    plt.legend(['四级四阶龙格-库塔法函数值数值计算结果',
                '四级四阶龙格-库塔法导数数值计算结果', '欧拉法函数值数值计算结果',
                '欧拉法导数数值计算结果'])
    plt.show()

import numpy as np
from matplotlib import pyplot as plt
plt.rcParams['font.sans-serif']=['SimHei']
plt.rcParams['axes.unicode_minus'] = False


def myfun1(x, y1, y2):
    return y2


def myfun2(x,y1, y2):
    return np.exp(2*x)*np.sin(x)-2*y1+2*y2


def ode(myfun1, myfun2, x, y1, y2, h, b):
    y1_r = []
    y2_r = []
    x_r = []
    y1_r.append(y1)
    y2_r.append(y2)
    x_r.append(x)
    for i in range(int(b / h)):
        k11 = h * myfun1(x, y1, y2)
        k12 = h * myfun2(x, y1, y2)
        k21 = h * myfun1(x + h / 2, y1 + 0.5 * k11, y2 + 0.5 * k12)
        k22 = h * myfun2(x + h / 2, y1 + 0.5 * k11, y2 + 0.5 * k12)
        k31 = h * myfun1(x + h / 2, y1 + 0.5 * k21, y2 + 0.5 * k22)
        k32 = h * myfun2(x + h / 2, y1 + 0.5 * k21, y2 + 0.5 * k22)
        k41 = h * myfun1(x + h, y1 +  k31, y2 + k32)
        k42 = h * myfun2(x + h, y1 +  k31, y2 + k32)
        y1 = y1 + (k11 + 2*k21 + 2*k31 + k41)/6
        y2 = y2 + (k12 + 2*k22 + 2*k32 + k42)/6
        y1_r.append(y1)
        y2_r.append(y2)
        x = x + h
        x_r.append(x)
    return x_r, y1_r, y2_r


x = np.linspace(0, 10, 20)
y = 0.2*np.exp(2*x)*(np.sin(x)-2*np.cos(x))
x_r, y1_r, y2_r = ode(myfun1, myfun2, 0, -0.4, -0.6, 0.01, 10)
plt.plot(x, y, 'ks:', markersize=3)
plt.plot(x_r, y1_r, 'b*:', markersize=1.5)
plt.plot(x_r, y2_r, 'g^:', markersize=1.5)
plt.legend(['函数值数值计算结果', '原函数', '导数数值计算结果'])
plt.show()

0到1区间

四级四阶龙格-库塔法和欧拉法的比较

import numpy as np
from matplotlib import pyplot as plt
plt.rcParams['font.sans-serif']=['SimHei']
plt.rcParams['axes.unicode_minus'] = False


def myfun1(x, y1, y2):
    return y2


def myfun2(x, y1, y2):
    return np.exp(2*x)*np.sin(x)-2*y1+2*y2


def ode1(myfun1, myfun2, x, y1, y2, h, b):
    y1_r = []
    y2_r = []
    x_r = []
    y1_r.append(y1)
    y2_r.append(y2)
    x_r.append(x)
    for i in range(int(b / h)):
        x_old = x
        y1_old = y1
        y2_old = y2
        y1 = y1 + h*myfun1(x_old, y1_old, y2_old)
        y2 = y2 + h*myfun2(x_old, y1_old, y2_old)
        y1_r.append(y1)
        y2_r.append(y2)
        x = x + h
        x_r.append(x)
    return x_r, y1_r, y2_r


def ode2(myfun1, myfun2, x, y1, y2, h, b):
    y1_r = []
    y2_r = []
    x_r = []
    y1_r.append(y1)
    y2_r.append(y2)
    x_r.append(x)
    for i in range(int(b / h)):
        k11 = h * myfun1(x, y1, y2)
        k12 = h * myfun2(x, y1, y2)
        k21 = h * myfun1(x + h / 2, y1 + 0.5 * k11, y2 + 0.5 * k12)
        k22 = h * myfun2(x + h / 2, y1 + 0.5 * k11, y2 + 0.5 * k12)
        k31 = h * myfun1(x + h / 2, y1 + 0.5 * k21, y2 + 0.5 * k22)
        k32 = h * myfun2(x + h / 2, y1 + 0.5 * k21, y2 + 0.5 * k22)
        k41 = h * myfun1(x + h, y1 + 0.5 * k31, y2 + k32)
        k42 = h * myfun2(x + h, y1 + 0.5 * k31, y2 + k32)
        y1 = y1 + (k11 + 2*k21 + 2*k31 + k41)/6
        y2 = y2 + (k12 + 2*k22 + 2*k32 + k42)/6
        y1_r.append(y1)
        y2_r.append(y2)
        x = x + h
        x_r.append(x)
    return x_r, y1_r, y2_r


x = np.linspace(0, 1, 20)
y = 0.2*np.exp(2*x)*(np.sin(x)-2*np.cos(x))
x_r, y1_r, y2_r = ode1(myfun1, myfun2, 0, -0.4, -0.6, 0.001, 1)
x_r1, y1_r1, y2_r1 = ode2(myfun1, myfun2, 0, -0.4, -0.6, 0.001, 1)
plt.plot(x, y, '1:', markersize=5)
plt.plot(x_r, y1_r, '2:', markersize=1.5)
plt.plot(x_r, y2_r, '3:', markersize=1.5)
plt.plot(x_r1, y1_r1, '4:', markersize=1.5)
plt.plot(x_r1, y2_r1, '8:', markersize=1.5)
plt.legend(['原函数', '四级四阶龙格-库塔法函数值数值计算结果',
            '四级四阶龙格-库塔法导数数值计算结果', '欧拉法函数值数值计算结果',
            '欧拉法导数数值计算结果'])
plt.show()