计算物理中matlab处理微分方程解析解和欧拉法数值解的算法演示

我先来看一个问题的引入：

我们根据题目给出的微分方程编写matlab求解代码如下：

syms cd m g v(t);
eqn = diff(v,t) == g - cd/m*v^2;
vt = dsolve(eqn,cond)

求解结果如下：

在得知相关初始条件后，对代码进一步设置求解：

syms cd m g v(t);                   % 定义符号变量
eqn = diff(v,t) == g - cd/m*v^2;    % 微分方程
cond = v(0) == 0;                   % 设置初始速度
vt = dsolve(eqn,cond);              % 解微分方程

vt = subs(subs(subs(vt,'m',68.1),'g',9.8),'cd',0.25);   % 迭代赋值

for i = 0:30
    vvt(i+1) = double(subs(vt,'t',i));  % 循环计算0~30s每隔1s的速度，并放入数组vvt
end

% 判断级数是否收敛
fn = vt;                       % 定义级数表达式（速度随时间变化的表达式，即微分方程的解）
sum = symsum(vt,t,0,inf);      % 对级数进行求和
vpa(sum)        % 输出近似值

%绘图
ft = 0:30;
plot(ft,vvt,'r-o','LineWidth',2)
title('前30s物体运动速度随时间的变化关系');
xlabel('时间t（s）');
ylabel('对应时刻的速度v（m/s）');
grid on; %打开栅格
hold on; %保持图像窗口，继续添加
plot(12,vvt(13),'b-+','LineWidth',2) %标明关注点
text(10,54,'12s时速度');

绘制的图像如下：

以上即为解析解图像，接下来我们利用欧拉法近似逼近得到数值解：
首先编写近似函数：

function vt = Numerical_solution(m,v0,cd,dt,t)
g=9.8;  % 重力加速度
for i = 0:dt:t
    v2 = v0 + (g-cd/m*v0^2)*dt; % v2为dt时间结束时的速度，v0为dt时间内开始时的速度
    vt(i+1,1) = i;              % 将循环0~ts内间隔dt的时刻保存在数组vt第一列
    vt(i+1,2) = v0;             % 将循环0~ts内间隔dt的速度保存在数组vt第二列
    v0 = v2;                    % 将上一次的末速度设置为下一次的初速度
end
plot(vt(:,1),vt(:,2),'b-','LineWidth',2);   % 绘图

与此同时在主脚本中调用函数，并将其与解析解绘制在同一图中：

syms cd m g v(t);                   % 定义符号变量
eqn = diff(v,t) == g - cd/m*v^2;    % 微分方程
cond = v(0) == 0;                   % 设置初始速度
vt = dsolve(eqn,cond);              % 解微分方程

vt = subs(subs(subs(vt,'m',68.1),'g',9.8),'cd',0.25);   % 迭代赋值

for i = 0:30
    vvt(i+1) = double(subs(vt,'t',i));  % 循环计算0~30s每隔1s的速度，并放入数组vvt
end

% 判断级数是否收敛
fn = vt;                       % 定义级数表达式（速度随时间变化的表达式，即微分方程的解）
sum = symsum(vt,t,0,inf);      % 对级数进行求和
vpa(sum)        % 输出近似值

%绘图
ft = 0:30;
plot(ft,vvt,'r-o','LineWidth',2)
title('前30s物体运动速度随时间的变化关系');
xlabel('时间t（s）');
ylabel('对应时刻的速度v（m/s）');
grid on; %打开栅格
hold on; %保持图像窗口，继续添加
plot(12,vvt(13),'b-+','LineWidth',2) %标明关注点
Numerical_solution(68.1,0,0.25,1,30); %数值解图像
text(10,54,'12s时速度');
legend('解析解','关注点','数值解','Location','best')

最终运行代码得到以下结果：

我们发现解析解与数值解存在一定的偏差，这是由于数值解采用的是近似逼近的求解算法，接下来我们将源代码脚本中的时间dt修改为0.1，即：
Numerical_solution(68.1,0,0.25,0.1,30); %数值解图像
与此同时避免函数文件报错，在此做如下修改：

function vt = Numerical_solution(m,v0,cd,dt,t)
g=9.8;  % 重力加速度
i = 0;
for di = 0:dt:t
    v2 = v0 + (g-cd/m*v0^2)*dt; % v2为dt时间结束时的速度，v0为dt时间内开始时的速度
    vt(i+1,1) = di;              % 将循环0~ts内间隔dt的时刻保存在数组vt第一列
    vt(i+1,2) = v0;             % 将循环0~ts内间隔dt的速度保存在数组vt第二列
    v0 = v2;                    % 将上一次的末速度设置为下一次的初速度
    i = i+1;
end
plot(vt(:,1),vt(:,2),'b-','LineWidth',2);   % 绘图  

重新运行代码得：

我们可以看到，随着数值解的求解精度提高，数值解的结果逐步趋近于解析解并最终与之重合，而在大部分实际应用问题中，很少存在解析解的情况，大多数情况都需要采取数值解近似逼近求解的方法！
而常用的数值解求解方法有：欧拉法公式，Runge-Kutta公式，Adams公式等，感兴趣的可以进一步去了解这些算法！