1.6 微分方程求解
很多工程问题，特别是力学类问题，常用微分方程来表示节点的位移偏移、应力等变化情况，这些微小的变化则采用数学上的微积分来表示，因此微分方程的应用求解，成为解决实际工程应用问题的有力解决方法，本节主要介绍数值积分运算以及微分方程数值解Ode等计算。
1.6.1 数值积分运算
数值积分是数值计算中常用的一个重要知识点，数值积分是一种在满足用户精度需求的情况下，对实际问题的近似处理。数值积分运算方法比较多，本小节主要讲述一种较常用的积分算法——龙贝格积分算法。
龙贝格积分算法计算步骤如下，其中
是一系列逼近原定积分的龙贝格积分值。
（1）计算首项
：
（2）对k＝1，2，3，…，计算如下：
对m＝1，2，…，k和i＝k，k-1，k=2，…1，计算：
上面的计算过程如表1-11所示。
表1-11 龙贝格积分计算表格
随着计算步骤的增加，
越来越逼近积分
。
在MATLAB中编程实现的龙贝格积分法的函数程序如下：
function [I,step]=Roberg(f,a,b,eps)
% 龙贝格积分法
% I: 积分值
% step: 积分划分的子区间次数
% f: 函数名
% a: 积分下限
% b: 积分上限
% eps: 积分精度
if(nargin==3) % 在缺省eps的情况下
eps=1.0e-4; % 默认值
end;
M=1; % 初始化
tol=10; % 前后两次积分差值
k=0; % 初始化
T=zeros(1,1); % 初始化
h=b-a; % 积分区间长度
% 计算首项
T(1,1)=(h/2)*(subs(sym(f),findsym(sym(f)),a)+subs(sym(f),findsym(sym(f)),b)); %初始值
while tol>eps
k=k+1; % 下标
h=h/2; % 区间长度
Q=0; % 初始化
for i=1:M
x=a+h*(2*i-1);
Q=Q+subs(sym(f),findsym(sym(f)),x); % 积分
end
T(k+1,1)=T(k,1)/2+h*Q; % 更新T值
M=2*M; % 2倍
for j=1:k
T(k+1,j+1)=T(k+1,j)+(T(k+1,j)-T(k,j))/(4^j-1); % 更新T值
end
tol=abs(T(k+1,j+1)-T(k,j)); % 前后两次积分差值
end
I=T(k+1,k+1);
step=k; % 积分划分的子区间次数
采用龙贝格积分法数值积分求解积分
和
。
在MATLAB命令窗口中输入下列命令：
>> [q,s]=Roberg('2*x^3',-2,3,1e-4)
q =
32.5000
s =
2
>> [q,s]=Roberg('2*x^1',-2,3,1e-4)
q =
5
s =
1
由龙贝格积分法可得到
。
1.6.2 微分方程数值解Ode
求解微分方程
的数值解，x取值区间[4,20]，初值条件
。
由
编写MATLAB函数文件如下：
function dy=fun(x,y)
dy=cos(x)+sin(y); % 函数
end
求解该函数的数值解，编程如下：
% Designed by Yu Shengwei From SWJTU University
% 2014年12月29日
clc,clear,close all % 清理命令区、清理工作区、关闭显示图形
warning off % 消除警告
feature jit off % 加速代码运行
tic % 运算计时
options = odeset('RelTol',1e-4,'AbsTol',1e-5); % Ode求解器参数设置
y0 = 1; % 初值
[T,Y] = ode45(@fun,[4 20],y0,options);
box on % 外框盒子边线
grid on % 网格化
axis equal % 轴相等
plot(T,Y,'linewidth',2) % 画图
toc % 计时结束'
运行程序输出图形如图1-39所示。
图1-39 数值解
求解微分方程的数值解常用的MATLAB函数调用如下：
说明如下。
（1）两个指令的调用格式相同，均为龙格库塔（Runge-Kutta）法。
（2）该指令是针对一阶常微分设计的。如果待解的方程是高阶微分方程，那么它必须先转化为形如
的一阶微分方程组，即“状态方程”。
（3）‘fun’是定义f（x，t）的函数名。该函数文件必须以
为一个列向量输出，以t，x为输入参量（t为“时间变量”，x为“状态变量”）。
（4）输入参量t0和tf分别是积分的起始值和终止值。
（5）输入参量x0为初始状态列向量。
（6）tol控制解的精度，可缺省。缺省时，ode23默认tol=1.e-3；ode45默认tol=1.e-6。
（7）输入参量trace控制求解的中间结果是否显示，可缺省，缺省时，默认为tol=0，表示不显示中间结果。
（8）输出参量t为“时间变量”，x为“状态变量”。
如求解著名的Van Der Pol方程
的数值解，并绘制其时间相应曲线和状态轨迹图。
首先令
，把
写成状态方程如下：
编写函数文件如下：
function xdot=VDP(t,x)
xdot=zeros(2,1); % 初始化，二元零向量
xdot(1)=(1-x(2)^2)*x(1)-x(2); % x1
xdot(2)=x(1); % x2
% xdot=[(1-x(2)^2)*x(1)-x(2);x(1)]; % 表示方法二
求解微分方程，编写主程序如下：
% Designed by Yu Shengwei From SWJTU University
% 2014年12月29日
clc,clear,close all % 清理命令区、清理工作区、关闭显示图形
warning off % 消除警告
feature jit off % 加速代码运行
tic % 运算计时
t0=0; % 时间初始值
tf=20; % 时间终止值
x0=[0,0.25]; % 初始值
[t,x]=ode45('VDP',[t0 tf],x0); % ode45龙格库塔求解
figure(1),
plot(t,x(:,1),':b',t,x(:,2),'-r') % 画图
legend('速度','位移') % 标记
figure(2), % 先建窗口画图
plot(x(:,1),x(:,2),'linewidth',2); % 画图
xlabel('速度');ylabel('位移')
toc % 计时结束'
运行程序输出结果如图1-40和图1-41所示。
图1-40 位移速度图形
图1-41 状态轨迹图
当然也可以将生成图1-40和图1-41的程序采用函数Function文件的形式给出，具体的程序如下：
function ysw1_17
% From SWJTU University
warning off % 消除警告
feature jit off % 加速代码运行
tic % 运算计时
t0=0; % 时间初始值
tf=20; % 时间终止值
x0=[0,0.25]; % 初始值
[t,x]=ode45(@VDP,[t0 tf],x0); % 求解
figure(1),
plot(t,x(:,1),':b',t,x(:,2),'-r') % 画图
legend('速度','位移')
figure(2), % 先建一个图形窗口
plot(x(:,1),x(:,2),'linewidth',2); % 画图
xlabel('速度');ylabel('位移')
toc % 计时结束'
end
function xdot=VDP(t,x)
xdot=zeros(2,1); % 初始化，二元零向量
xdot(1)=(1-x(2)^2)*x(1)-x(2); % x1
xdot(2)=x(1); % x2
% xdot=[(1-x(2)^2)*x(1)-x(2);x(1)]; % 表示方法二
end