我先来看一个问题的引入：
我们根据题目给出的微分方程编写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:30vvt(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:tv2 = 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:30vvt(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:tv2 = 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公式等，感兴趣的可以进一步去了解这些算法！
