Matlab求解微分方程（解析解与数值解）

Matlab求微分方程解析解

dsolve(eqns,conds,options) eqns:微分方程（组）、conds：初值条件、options:精度设置

例题1

%例题一：y-Dy=2x(Dny表示y的n阶导),y(0)=3
dsolve('y-Dy=2*x','y(0)=3','x')

这种求解形式虽然简单，但是警告内容说它会在未来版本中被移除，因此推荐使用下面的方法使用dsolve函数

syms y(x)
eqn1=(y-diff(y,x)==2*x)
cond1=(y(0)==3)
dsolve(eqn1,cond1)

例题2

%例题二：D2y+4Dy+29y=0,y(0)=0,Dy(0)=15
syms y(x)
eqn2=(diff(y,x,2)+4*diff(y,x)+29*y==0)
Dy=diff(y,x)
cond2=[(y(0)==0),(Dy(0)==15)]
dsolve(eqn2,cond2)

例题3

syms x(t) y(t) z(t)
eqn3=[(diff(x,t)==2*x-3*y+3*z+t),(diff(y, t)==4*x-5*y+3*z+t),(diff(z,t)==4*x-4*y+2*z+t)]
ans3=dsolve(eqn3)

Matlab求微分方程数值解

一阶微分方程

[x,y]=solver('f',ts,x0,options)

①solver：求解器（以具体求解器作为函数名）
最常用的两个求解器包括ode45和ode15s。ode45采用4-5阶Runge-Kutta法，适用于求解非刚性问题（求解出的函数图像不发生突变）；ode15s基于数值差分公式，适用于求解刚性问题（求解出的函数图像发生突变）。
因此在实际操作中，通常先采用ode45作为求解器，如果运行过程卡顿，说明面临刚性问题，应切换ode15s。

②f:包含微分方程的函数文件名f.m，除了使用字符串‘f’外可传入@f，标准形式dy=f（x,y）(dy;一阶微分)

③ts=[$t_a,t_b$]，表示自变量的范围，可使用$t_a:step:t_b$设置步距

④x0:函数的初值，n个函数对应n个初值

⑤options:设置求解精度，options=odeset(‘reltol’,rt,‘abstol’,at),默认相对误差rt=1e-3，绝对误差at=1e-6

⑥[x,y] x:返回所取自变量的离散值 y：返回各个函数对应于x的函数值

例题一

%df1.m
function dy=df1(x,y) 
    dy=y-2*x;
end

%练习一
[x1,y1]=ode45(@df1,[0,2],3)
figure(1)
plot(x1,y1,'r*-')

例题二

function dy=df2(x,y)
   dy=zeros(3,1);
   dy(1)=y(2)*y(3);
   dy(2)=-y(1)*y(3);
   dy(3)=-0.51*y(1)*y(2);
end

%练习二
[x2,y2]=ode45(@df2,[0 4*pi],[0 1 1])
figure(2)
plot(x2,y2(:,1),'o',x2,y2(:,2),'*',x2,y2(:,3),'+')
legend('y1','y2','y3')
axis([0 4*pi -inf inf])

高阶微分方程

处理高阶微分方程，需要先通过变量转化将问题转变为一阶微分方程，然后使用solver求解

例题

$设置y_1=y^{{}^{\prime }},y_2=y,则（1+x^2）y_1^{{}^{\prime }}=2x{y}_1,y_2^{{}^{\prime }}=y_1,二阶微分方程转化为一阶微分方程组$

%df3.m
function dy=df3(x,y)
    dy=zeros(2,1)
    dy(2)=y(1)
    dy(1)=(2*x*y(1))/(1+x^2)
end

%高阶微分方程例题 y1=y',y2=y
[x,y]=ode45(@df3,[-2,2],[4 3])
figure(1)
plot(x,y(:,2),'r*')