编程手之数学建模学习——微分方程

分为符号解和数值解。

符号解是精确的，可以表示为表达式或函数；

数值解使用离散的数值数据来近似微分方程的解；

一、符号解求解

微分方程求解符号解：

采用dsolve函数：

clc, clear, syms y(x)  %定义符号函数ｙ，自变量为x
dy=diff(y);  %定义y的一阶导数，目的是下面赋初值
y=dsolve(diff(y,2)-2*diff(y)+y==exp(x), y(0)==1, dy(0)==-1)

微分方程组求解符号解：

矩阵形式（一阶微分方程组）

clc,clear
syms x(t)  [3,1]  %定义符号向量函数，x(t)后面有空格
A=[1,0,0;2,1,-2;3,2,1]; B=[0;0;exp(t)*cos(2*t)];
x0=[0;1;1]; %初值条件
[s1,s2,s3]=dsolve(diff(x)==A*x+B,x(0)==x0) %求符号解

 组合形式（高阶微分方程组）

clc,clear, syms f(x) g(x)  %定义符号函数
df=diff(f);  %定义f的一阶导数
[s1,s2]=dsolve(diff(f,2)+g-2==1, diff(g)+diff(f)==1,...
        df(1)==0, f(0)==0, g(0)==0)

二、数值解求解（重点）

初值问题求解数值解：

 ode45的调用 （必须转为一次方程组）

clc, clear, close all
dy=@(x,y)[y(2); sqrt(1+y(2)^2)/5/(1-x)];
[x,y]=ode45(dy,[0,0.999999],[0;0])
plot(x, y(:,1)), xlabel('$x$','Interpreter','Latex')
ylabel('$y$','Interpreter','Latex','Rotation',0)

 deval函数确定特定点上的值

这里的s也就是前面ode45得到的y和s。

 边值问题求解数值解：

给出自变量区域边界的函数值，求解微分方程的数值解。

采用bvpinit函数给出初始猜测解具体数值，采用bvp4c计算数值解。

bvpinit(自变量取值，猜测符号解)；bvp4c(一阶微分方程组，边值条件，猜测初始解具体数值)

clc, clear
eq=@(x,y,mu)[y(2);-mu*y(1)]; %定义一阶方程组的匿名函数
bd=@(ya,yb,mu)[ya(1);ya(2)-1;yb(1)+yb(2)];  %定义边值条件的匿名函数
guess=@(x)[sin(2*x);2*cos(2*x)]; %定义初始猜测解的匿名函数
guess_structure=bvpinit(linspace(0,1,10),guess,5); %给出初始猜测解的结构,mu=5
sol=bvp4c(eq,bd,guess_structure); %计算数值解
plot(sol.x, sol.y(1,:), '-', sol.x, sol.yp(1,:), '--', 'LineWidth',2)
xlabel('$x$','FontSize',12,'Interpreter','Latex')
ylabel('$y$','FontSize',12,'Interpreter','Latex','Rotation',0)
legend('$y_1$','$y_2$','Interpreter','Latex')