数据拟合 龙格库塔 数学外推

Question

要求使用四阶龙格库塔法求解四个方程的系数a，b和C，并作出四个y的数据拟合曲线（分开画）

Answer 1

以下答案由GPT-3.5大模型与博主波罗歌共同编写：
首先，关于四阶龙格库塔法的求解方法，请参考链接：https://zh.wikipedia.org/wiki/%E9%BE%99%E6%A0%BC%EF%BC%8D%E5%BA%93%E5%A1%94%E6%B3%95

针对题目中的需求，我们需要对四个方程进行求解，并作出数据拟合曲线。具体步骤如下：

1. 定义四个方程，根据实际问题确定。

2. 接下来需要使用四阶龙格库塔法求解方程。具体实现方法可以参考以下代码：

import numpy as np

# 定义四个方程，这里只是示例，请根据实际问题进行修改
def f1(x, y1, y2, y3, y4):
    return 2 * y2 + 3 * y3 - 6 * y4

def f2(x, y1, y2, y3, y4):
    return -2 * y1 + y2 - 2 * y3 + y4 + np.sin(x)

def f3(x, y1, y2, y3, y4):
    return 7 * y1 - 8 * y2 + 9 * y3 - 10 * y4 + x**2

def f4(x, y1, y2, y3, y4):
    return -y1 + 3 * y2 - 4 * y3 + 2 * y4 + np.exp(x)

# 定义求解过程
def runge_kutta4(f1, f2, f3, f4, a, b, h, y1, y2, y3, y4):
    n = int((b - a) / h) # 计算步数
    X = [a] # 存储X的值
    Y1 = [y1] # 存储Y1的值
    Y2 = [y2] # 存储Y2的值
    Y3 = [y3] # 存储Y3的值
    Y4 = [y4] # 存储Y4的值
    for i in range(n):
        k11 = h * f1(X[-1], Y1[-1], Y2[-1], Y3[-1], Y4[-1])
        k12 = h * f2(X[-1], Y1[-1], Y2[-1], Y3[-1], Y4[-1])
        k13 = h * f3(X[-1], Y1[-1], Y2[-1], Y3[-1], Y4[-1])
        k14 = h * f4(X[-1], Y1[-1], Y2[-1], Y3[-1], Y4[-1])
        
        k21 = h * f1(X[-1] + h/2, Y1[-1] + k11/2, Y2[-1] + k12/2, Y3[-1] + k13/2, Y4[-1] + k14/2)
        k22 = h * f2(X[-1] + h/2, Y1[-1] + k11/2, Y2[-1] + k12/2, Y3[-1] + k13/2, Y4[-1] + k14/2)
        k23 = h * f3(X[-1] + h/2, Y1[-1] + k11/2, Y2[-1] + k12/2, Y3[-1] + k13/2, Y4[-1] + k14/2)
        k24 = h * f4(X[-1] + h/2, Y1[-1] + k11/2, Y2[-1] + k12/2, Y3[-1] + k13/2, Y4[-1] + k14/2)
        
        k31 = h * f1(X[-1] + h/2, Y1[-1] + k21/2, Y2[-1] + k22/2, Y3[-1] + k23/2, Y4[-1] + k24/2)
        k32 = h * f2(X[-1] + h/2, Y1[-1] + k21/2, Y2[-1] + k22/2, Y3[-1] + k23/2, Y4[-1] + k24/2)
        k33 = h * f3(X[-1] + h/2, Y1[-1] + k21/2, Y2[-1] + k22/2, Y3[-1] + k23/2, Y4[-1] + k24/2)
        k34 = h * f4(X[-1] + h/2, Y1[-1] + k21/2, Y2[-1] + k22/2, Y3[-1] + k23/2, Y4[-1] + k24/2)
        
        k41 = h * f1(X[-1] + h, Y1[-1] + k31, Y2[-1] + k32, Y3[-1] + k33, Y4[-1] + k34)
        k42 = h * f2(X[-1] + h, Y1[-1] + k31, Y2[-1] + k32, Y3[-1] + k33, Y4[-1] + k34)
        k43 = h * f3(X[-1] + h, Y1[-1] + k31, Y2[-1] + k32, Y3[-1] + k33, Y4[-1] + k34)
        k44 = h * f4(X[-1] + h, Y1[-1] + k31, Y2[-1] + k32, Y3[-1] + k33, Y4[-1] + k34)
        
        Y1.append(Y1[-1] + 1/6 * (k11 + 2*k21 + 2*k31 + k41))
        Y2.append(Y2[-1] + 1/6 * (k12 + 2*k22 + 2*k32 + k42))
        Y3.append(Y3[-1] + 1/6 * (k13 + 2*k23 + 2*k33 + k43))
        Y4.append(Y4[-1] + 1/6 * (k14 + 2*k24 + 2*k34 + k44))
        X.append(X[-1] + h)
        
    return X, Y1, Y2, Y3, Y4

1. 定义要拟合的数据点集合，这里假设我们有一组数据点(x,y)，将其存储到两个列表中，分别表示x和y。

2. 使用上面定义的runge_kutta4函数求解对应的四个方程，得到四组y值，存储到一个二维数组matrix中。

3. 对每组数据进行拟合，得到对应的拟合曲线。可以使用numpy.polyfit函数进行多项式拟合。下面是示例代码：

import matplotlib.pyplot as plt

# 定义拟合阶数
n = 5

# 准备数据点
x = [0.0, 0.1, 0.2, 0.3, 0.4, 0.5]
y = [0.0, 0.3, 0.55, 0.6, 0.45, 0.18]

# 求解四个方程，得到四组y值
X, Y1, Y2, Y3, Y4 = runge_kutta4(f1, f2, f3, f4, 0.0, 0.5, 0.1, 0.0, 0.0, 0.0, 0.0)
matrix = np.array([Y1, Y2, Y3, Y4]).T

# 绘制拟合曲线
fig, axs = plt.subplots(nrows=2, ncols=2)
axs[0, 0].plot(x, y, 'o')
fit = np.polyfit(x, y, n)
equation = np.poly1d(fit)
x_axis = np.linspace(0, 0.5, 100)
axs[0, 0].plot(x_axis, equation(x_axis), '-')
axs[0, 0].set_title('y1')

axs[0, 1].plot(x, y, 'o')
fit = np.polyfit(x, matrix[:, 0], n)
equation = np.poly1d(fit)
axs[0, 1].plot(x_axis, equation(x_axis), '-')
axs[0, 1].set_title('y2')

axs[1, 0].plot(x, y, 'o')
fit = np.polyfit(x, matrix[:, 1], n)
equation = np.poly1d(fit)
axs[1, 0].plot(x_axis, equation(x_axis), '-')
axs[1, 0].set_title('y3')

axs[1, 1].plot(x, y, 'o')
fit = np.polyfit(x, matrix[:, 2], n)
equation = np.poly1d(fit)
axs[1, 1].plot(x_axis, equation(x_axis), '-')
axs[1, 1].set_title('y4')

plt.show()

上面的代码中，我们首先绘制了原始数据点的图形，并对每组数据进行了拟合，得到对应的拟合曲线。其中，np.polyfit函数的第一个参数是x值，第二个参数是y值，第三个参数是拟合阶数。最后使用np.poly1d函数创建一个一维多项式函数，参数为拟合系数。
如果我的回答解决了您的问题，请采纳！