【数学建模】14 微分方程模型求解方法

1 MATLAB数值微积分

（1）差分与微分 • taylor 符号泰勒展开 • polyder 多项式求导 • diff 数值差分或符号求导

dx = diff(x) %返回向量x的差分

• gradient 数值梯度

Fx = gradient(F,x) % 返回向量F表示的一元函数沿x方向的数值梯度（F'(x)）,其中，x是与同维度的向量.x是自变量，F是因变量
[Fx,Fy] = gradient(F,x,y)%返回矩阵F表示的二元函数的数值梯度(F'x,F'y),当F为m*n矩阵时x、y分别为n维和m维的向量。x和y是因变量

举例，求梯度，两种方式 第一种：通过x和y的差分比.误差比较大 第二种：用gradient函数求。两端的数据偏大，中间的数据比较准确

clear
x = [1 1.1 1.2 1.3];
y = x.^3;
%第一种方式
dy = diff(y)./diff(x) 
%第二种方式
dy = gradient(y,x)

（2）积分 • trapz：梯形法积分

Z = trapz(x,y)% x是表示积分区间的离散化向量，y是与x同维数的向量，表示被积函数；z返回积分的近似值。

• dblquad：矩形区域的二重积分(二重积分)

Z = dblquad('Fun',a,b,c,d)
%其中：Fun表示被积函数的M函数名，a、b表示x的上下限，c,d表示x的上下限

• quad：变步长数值积分

Z = quad('Fun',A,B,Tol)
% 其中：Fun表示被积函数的M函数名，A、B上下限，Tol精度，缺省值为1e-3

• int：符号积分

int(s) %符号表达式s的不定积分
int(s,v)% 符号表达式s关于变量v的不定积分
int(s,a,b)%表达式s的定积分，a，b分别为下限和上限

• quad8：高精度数值积分 举例1，计算积分： z = ∫ − 1 1 e − x 2 d x z = \int^{1}_{-1}{e^{-x^2}dx} z=∫−11​e−x2dx

clear
x = -1:0.1:1;
y = exp(-x.^2);
% 方法一
trapz(x,y)
% 方法二
quad(@(x)exp(-x.^2),-1,1)

举例2，计算二重积分 z = ∫ 0 1 d x ∫ 0 5 s i n ( x y ) d y z = \int^{1}_{0}{dx}\int^{5}_{0}{sin(xy)dy} z=∫01​dx∫05​sin(xy)dy

Z = dblquad(@(x,y)sin(x.*y),0,1,0,5)

举例3，计算积分 z = ∫ 0 1 e − x + s i n ( x ) d x z = \int^{1}_{0}{e^{-x} + sin(x) dx} z=∫01​e−x+sin(x)dx

t2 = int(exp(-x) + sin(x),0,1)
t2 = vpa(t2) %求近似解

2 微分方程数值解

3 MATLAB求解常微分方程

• ode23 二、三阶Runge-Kutta法 • ode15s 刚性方程组解法 • ode45 四、五阶Runge-Kutta法 • dsolve 符号解析解

[tout,yout] = ode45('yprime',tspan,y0)
%这是字符串yprime是用以表示f(t,y)的M文件名
%tspan = [t0,tf] 表示自变量初始值t0和终值tf.y0表示初值向量y0。
%输出列向量tout表示节点（t0,t1,..,tn）
%输出矩阵yout表示数值解，每一列对应y的一个分量。

举例1，解方程 y ′ = t − 2 t y y ( 0 ) = 1 , 0 < t < 1 y' = t- \frac{2t}{y} \\y(0) = 1,0