新手求助：如何用Matlab ODE23求解牛顿冷却微分方程

嘿，作为曾经也对着MATLAB ODE函数摸不着头脑的过来人，我来一步步带你搞定这个牛顿冷却方程的求解！咱们先把问题再理清楚：你要解的微分方程是：

dT_c/dt = -r(T_c - T_s)
其中参数：T_s=19（环境温度），初始温度T_c(0)=84，r=0.025，求解时间区间是[0, 300]秒。

步骤1：理解ode23的输入要求 #

MATLAB的ode23是用来求解一阶常微分方程（组）的函数，它的基本调用格式是：

[t, y] = ode23(ode_function, tspan, y0);

• ode_function：你要解的微分方程对应的函数，必须接受两个输入参数：时间t和当前状态变量y，输出是状态变量的导数dy/dt
• tspan：求解的时间区间，格式为[起始时间 结束时间]
• y0：状态变量的初始值

你的方程里，状态变量就是T_c，所以y=T_c，dy/dt就是dT_c/dt。

步骤2：定义你的ODE函数（两种可选方式） #

方式一：匿名函数（适合简单方程，无需单独写文件） #

因为你的方程里所有参数都是常数，直接在命令行或者脚本里定义匿名函数就行：

% 定义微分方程：dTc/dt = -0.025*(Tc - 19)
dTc_dt = @(t, Tc) -0.025*(Tc - 19);

这里t虽然在方程里没直接用到，但ode23要求必须把它作为第一个参数传入，所以不能省略。

方式二：单独的函数文件（适合复杂方程，方便复用） #

如果你以后要修改参数或者扩展方程，写一个单独的.m文件更方便。比如创建名为newton_cooling.m的文件，内容如下：

function dTc_dt = newton_cooling(t, Tc)
    % 定义方程参数
    T_s = 19;
    r = 0.025;
    % 计算温度的导数
    dTc_dt = -r*(Tc - T_s);
end

同样，t作为第一个参数必须保留，即使方程里没用到它。

步骤3：设置初始条件和时间区间 #

在脚本或者命令行里定义初始温度和时间区间：

Tc0 = 84;          % 初始温度
tspan = [0, 300];  % 求解的时间范围：0到300秒

步骤4：调用ode23求解 #

根据你选的函数方式，调用ode23：

• 如果用匿名函数：

[t, Tc] = ode23(dTc_dt, tspan, Tc0);

• 如果用单独函数文件：

[t, Tc] = ode23(@newton_cooling, tspan, Tc0);

这里的输出：

• t：ode23自动选取的时间点数组（自适应步长）
• Tc：对应每个时间点的温度数组

步骤5：可视化求解结果（推荐） #

为了直观看到温度变化，咱们画个图：

plot(t, Tc, 'b-', 'LineWidth', 2);
xlabel('时间 (秒)');
ylabel('温度 (℃)');
title('牛顿冷却过程温度变化');
grid on;
hold on;
% 画环境温度的参考线
plot(tspan, ones(size(tspan))*19, 'r--', 'LineWidth', 1.5);
legend('物体温度', '环境温度');
hold off;

运行后你会看到温度逐渐趋近于环境温度19℃，符合牛顿冷却定律的预期。

可选：验证数值解的正确性 #

牛顿冷却方程有解析解：

Tc(t) = T_s + (Tc0 - T_s)exp(-rt)
你可以把数值解和解析解对比，看看误差：

% 计算解析解
Tc_analytical = 19 + (84 - 19)*exp(-0.025*t);
% 绘制误差曲线
figure;
plot(t, abs(Tc - Tc_analytical), 'g-', 'LineWidth', 1.5);
xlabel('时间 (秒)');
ylabel('数值解与解析解的误差');
title('ODE23求解误差');
grid on;

你会看到误差非常小，说明ode23的求解结果是可靠的。

内容的提问来源于stack exchange，提问作者user3440639