数学建模——微分方程

目录

常微分方程符号解

常微分方程的数值解

初值问题

边值问题

 其他

---

常微分方程符号解

 diff(函数,n)%求函数的n阶导数

dsolve(方程1,方程2,...,方程n,初始条件,自变量)

simplify(s) 对表达式s进行化简

clc,clear
syms y(x);
dy = diff(y);
s1 = dsolve((1-x)*diff(y,2)==sqrt(1+dy^2)/5,y(0)==0,dy(0)==0,x)
simplify(s1)

clc,clear
syms y(x);
out = dsolve(diff(y)==-2*y+2*x^2+2*x,y(0)==1)

clc,clear
syms y(t);%同时声明了y和t
dy=diff(y);d2y=diff(y,2);d3y=diff(y,3);%定义导数是为了后面幅值
u =exp(-t)*cos(t);
y = dsolve(diff(y,4)+10*diff(y,3)+35*diff(y,2)+50*diff(y)+...
    24*y==diff(u,2),y(0)==0,dy(0)==-1,d2y(0)==1,d3y(0)==1)

clc,clear
syms x(t) [3,1];%定义符号向量函数，x(t)后面要有空格
A = [3,-1,1;2,0,-1;1,-1,2];%定义系数矩阵
[s1,s2,s3]=dsolve(diff(x)==A*x,x(0)==[1;1;1])

常微分方程的数值解

初值问题

一阶微分方程初值问题的一般格式：

常用ode45进行求解，有两种调用格式

 fun是匿名函数，tspan=[t0,tfinal]是求解区间，y0是初值，t、y是所求区间的对应点和该点的函数值，s是一个结构数组。tspan的t0必须是初始条件中自变量的取值，tfinal可以比t0小。利用函数deval可以计算区间内任意点的函数值：

y = deval(s,x);

 例子：

clc,clear,close all
syms y(x);
%求符号解
y = dsolve(diff(y)==-2*y+2*x^2+2*x,y(0)==1)
%求数值解
dy = @(x,y)-2*y+2*x^2+2*x;%匿名函数第一个应该是自变量，第二个是因变量
[sx,sy]=ode45(dy,[0,0.5],1);    
fplot(y,[0,0.5])%画符号解的图像
hold on 
plot(sx,sy,'*')%画数值解的图像
legend('符号解','数值解')

matlab不能求高阶微分方程，必须化为一阶微分方程组

例子：

应该先化为一阶微分方程组：

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))

clc,clear,close all
m = 2;
df = @(x,f,m)[f(2);m*f(2)-f(1)+exp(x)];%注意是列向量
s1 = ode45(@(x,f)df(x,f,m),[0,1],[1;-1]);
s2 = ode45(@(x,f)df(x,f,m),[0,-1],[1;-1]);
fplot(@(x)deval(s1,x,1),[0,1],'-ok');hold on
fplot(@(x)deval(s2,x,1),[-1,0],'-*k');

 只能求取一阶，所以定义两个未知数，ode45参数初值必须是第二个参数左端值对应的函数值，初值对应横轴0，所以将坐标范围分解为[0,-1]，[0,1]。匿名函数可以有三个参数，但ode45函数里需重新定义匿名函数，只能有两个形参。

边值问题

初值问题指定解条件的自变量是同一个点

 边值问题指定解条件的自变量是多个点

 该问题的微分方程标准形式：

 

 bc指边界条件，边界条件中自变量的取值也是求解的区间，即求解区间为[a,b]，p值有关的参数

相关函数：

sol = bvp4c(odefun,bcfun,solinit,options,p1,p2,…)

odefun 函数为微分方程的函数句柄。bcfun 函数为微分方程的边界条件句柄，将边界条件所有项移到左边，右边为0。solinit为包含初始估计解的结构体，应使用bvpinit函数创建，solinit=bvpinit(x,y)，x向量为初始网格的排序节点，应满足边界条件a= solinit.x(1)且b =solinit.x(end)，向量y为初始估计解，solinit.y(:,i)为在solinit.x(i)节点处的初始估计解。输出参数sol 为包含数值解的一个结构：

 为使曲线变得更光滑，需要在中间插入一些点，可以使用deval函数。

例子：

 化为一阶微分方程组：

 作为初始猜测解，可以随便取

clc,clear,close all
solinit = bvpinit(linspace(-1,1,20),@dropinit);
sol = bvp4c(@drop,@dropbc,solinit)
fill(sol.x,sol.y(1,:),[0.7 0.7 0.7]);
axis([-1 1 0 1]);

function yprime = drop(x,y) %微分方程组
yprime = [y(2),(y(1)-1)*(1+y(2)^2)^(3/2)];end

function res = dropbc(ya,yb)%边界条件
res = [ya(1);yb(1)];end     %y(1)表示y的0阶导，y(2)表示一阶导数，ya表示左边界，意思是ya(1)=0

function yinit = dropinit(x)%猜测函数
yinit = [x.^2,2*x];end

 例子：

 化为一阶微分方程组：

 

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);
sol = bvp4c(eq,bd,guess_structure);
plot(sol.x,sol.y(1,:),'-',sol.x,sol.yp(1,:),'--')

 其他

fminbnd用于查找单变量函数在定区间上的最小值