MATLAB实验-数学建模

编写程序，利用4阶Runge-Kutta方法求y’-y+2t/y=0,y(0)=1,0<t<4的数值解，步长h取0.1，在同一图形窗口画出数值解的曲线和精确解的曲线，添加标题、图例；计算在各个节点处的平均绝对误差、最大绝对误差。

MATLAB 中实现 4 阶 Runge-Kutta 方法的基本步骤

定义微分方程

首先，定义要求解的微分方程。在 MATLAB 中，创建一个函数，该函数接受当前时间点 t 和当前函数值 y 作为输入，并返回微分方程的导数值。

function dydt = odefun(t, y)
    % 微分方程的定义
    dydt = y - 2*t./y;
end

实现 4 阶 Runge-Kutta 方法

接着，实现 4 阶 Runge-Kutta 方法来求解微分方程。该方法在每个步骤中使用当前和预测的斜率值来估计函数的下一个值。

function [t, y] = runge_kutta(odefun, y0, tspan)
    % 4 阶 Runge-Kutta 方法的实现
    % 这里 odefun 是微分方程函数，y0 是初始值，tspan 是时间范围
end

完整matlab代码：

function runge_kutta_example
    % 初始条件
    y0 = 1;
    
    % 时间范围
    tspan = 0:0.1:4;
    
    % 使用 4 阶 Runge-Kutta 方法求解
    [t, y_rk4] = runge_kutta(@odefun, y0, tspan);
    
    % 精确解
    y_exact = sqrt(2*t + 1);

    % 计算平均绝对误差和最大绝对误差
    average_absolute_error = mean(abs(y_rk4 - y_exact));
    maximum_absolute_error = max(abs(y_rk4 - y_exact));
    
    % 绘制结果
    figure;
    plot(t, y_rk4, 'b-o', t, y_exact, 'r-');
    title('Runge-Kutta 4 阶近似解 vs 精确解');
    xlabel('t');
    ylabel('y');
    legend('RK4 近似解', '精确解');
    grid on;

    % 显示误差
    disp(['平均绝对误差: ', num2str(average_absolute_error)]);
    disp(['最大绝对误差: ', num2str(maximum_absolute_error)]);
end

function dydt = odefun(t, y)
    % 微分方程函数
    dydt = y - 2*t./y;
end

function [t, y] = runge_kutta(odefun, y0, tspan)
    % 4 阶 Runge-Kutta 方法实现
    n = length(tspan);
    y = zeros(n, 1);
    y(1) = y0;
    h = tspan(2) - tspan(1);

    for i = 1:(n-1)
        k1 = odefun(tspan(i), y(i));
        k2 = odefun(tspan(i) + 0.5*h, y(i) + 0.5*h*k1);
        k3 = odefun(tspan(i) + 0.5*h, y(i) + 0.5*h*k2);
        k4 = odefun(tspan(i) + h, y(i) + h*k3);
        y(i+1) = y(i) + (h/6)*(k1 + 2*k2 + 2*k3 + k4);
    end

    t = tspan;
end

运行结果：

 

 

编写程序，利用改进欧拉法求y’-y+2t/y=0,y(0)=1,0<t<4的数值解，步长h分别取0.1，0.01，在同一图形窗口中画出2条数值解的曲线和精确解的曲线，添加标题、图例；计算在各个节点处的平均绝对误差、最大绝对误差。提示：实际计算时，采用预估-校正公式进行计算。

要在 MATLAB 中使用改进欧拉法（有时也称为启发式欧拉法或修正欧拉法）求解给定的常微分方程，基本的步骤：

1. 定义微分方程 y′−y+2t​/y=0。
2. 编写一个函数实现改进欧拉法。
3. 对于给定的步长（h = 0.1 和 h = 0.01），使用该方法计算数值解。
4. 计算精确解 y(t)=sqrt(1+2t^2)。
5. 在同一图形窗口中绘制数值解和精确解的曲线。
6. 计算平均绝对误差和最大绝对误差。

 MATLAB 代码：

function [t, y] = improved_euler_method(dydt, t_range, y0, h)
    t = t_range(1):h:t_range(2);
    y = zeros(1, length(t));
    y(1) = y0;
    for i = 1:(length(t) - 1)
        f1 = dydt(t(i), y(i));
        y_pred = y(i) + f1 * h;
        f2 = dydt(t(i+1), y_pred);
        y(i+1) = y(i) + (h/2) * (f1 + f2);
    end
end

% 微分方程的右侧函数
function dy = dydt(t, y)
    dy = y - 2 * t / y;
end

% 给定的参数
t_range = [0, 4];
y0 = 1;
h1 = 0.1;
h2 = 0.01;

% 使用改进欧拉法计算数值解
[t1, y1] = improved_euler_method(@dydt, t_range, y0, h1);
[t2, y2] = improved_euler_method(@dydt, t_range, y0, h2);

% 计算精确解
t_exact = linspace(t_range(1), t_range(2), 1000);
y_exact = sqrt(1 + 2 * t_exact.^2);

% 绘图
figure;
plot(t1, y1, 'b.-', 'DisplayName', ['Improved Euler h=' num2str(h1)]);
hold on;
plot(t2, y2, 'g.-', 'DisplayName', ['Improved Euler h=' num2str(h2)]);
plot(t_exact, y_exact, 'r', 'DisplayName', 'Exact Solution');
title('Comparison of Improved Euler Method and Exact Solution');
xlabel('t');
ylabel('y(t)');
legend;
grid on;

% 计算误差
y_exact1 = sqrt(1 + 2 * t1.^2);
y_exact2 = sqrt(1 + 2 * t2.^2);
errors1 = abs(y1 - y_exact1);
errors2 = abs(y2 - y_exact2);

avg_abs_error1 = mean(errors1);
max_abs_error1 = max(errors1);
avg_abs_error2 = mean(errors2);
max_abs_error2 = max(errors2);

fprintf('For h=%0.2f, Average Absolute Error: %0.5f, Max Absolute Error: %0.5f\n', h1, avg_abs_error1, max_abs_error1);
fprintf('For h=%0.2f, Average Absolute Error: %0.5f, Max Absolute Error: %0.5f\n', h2, avg_abs_error2, max_abs_error2);

运行结果图：