MATLAB微分方程建模------2019/7/22

最新推荐文章于 2023-11-29 23:14:15 发布

立里∑

最新推荐文章于 2023-11-29 23:14:15 发布

阅读量2.8k

点赞数 3

分类专栏：微分方程建模

本文链接：https://blog.csdn.net/qq_41218103/article/details/96825125

版权

微分方程建模专栏收录该内容

2 篇文章 0 订阅

订阅专栏

微分方程建模

将实际问题转化为微分问题的步骤：
（1）根据实际要求确定要研究的量（自变量，未知函数，必要的参数等）并确定坐标系
（2）找出这些量所满足的基本规律
（3）运用规律列出方程和定解条件：
1）按规律直接列方程
2）微元分析法与任意区域上取积分的方法
3）模拟近似法
MATLAB求微分方程的符号解
（1）命令：
$[y1,y2,\cdots,yn]=dsolve(eqns,conds,name,value)$
其中：eqns为符号微分方程（组）；conds为初值条件或边值条件；name，value为可选的成对参数。

（2）求解线性常微分方程组
1）一阶齐次线性微分方程组
$X'=AX,X=\begin{bmatrix}x_{1}\\ \vdots\\ x_{n}\end{bmatrix},A=\begin{bmatrix}a_{11} \cdots a_{1n} \\ \vdots \ddots \vdots\\a_{n1}\cdots a_{nn} \end{bmatrix}$

式中： $^{'}$ 为对 $t$ 求导数， $e^{At}$ 是它的基解矩阵，解为 $X(t)=e^{A(t-t_{0})}X_{0}$
$d s o l v e$ 可以直接求解齐次线性微分方程组 $X^{'} = A X$
2)非齐次线性方程组
由参数变易法可求得初值问题
$X'=AX+f(t),X(t_{0})=X_{0}$
的解为：
$X(t)=e^{A(t-t_{0})}X_{0}+\int_{t_{0}}^{t} e^{A(t-s)f(s)}f(s)ds$
3. 初值问题的MATLAB数值解

简单的一阶方程的初值问题：
$\left\{ \begin{aligned} y'=f(x,y) \\ y(x_{0})=y_{0} \end{aligned} \right.$
MATLAB的函数形式是： $[t, y] = s o l v e r (^{'} F^{'}, t s p a n, y 0)$ ：
其中：
$F$ 是 $y^{'} = f (x, y)$ 的右端函数, $t s p a n = [t 0, t f i n a l]$ 是求解区间， $y 0$ 是初值。
$s o l v e r$ 为 $o d e 45, o d e 23, o d e 113$ 中的一个

检验方法：可以在代码中先用ode23求解的结果与用ode45求的结果进行对比，若差不多表明求解正确。

例6.7

yx=@(x,y) -2*y+2*x^2+2*x;  %定义微分方程右端的匿名函数
[x,y]=ode45(yx,[0,0.5],1)  %第一种返回格式
so1=ode45(yx,[0,0.5],1)    %第二种返回格式
y2=deval(so1,x)            %利用deval命令和返回的结构数组so1来计算x点的函数值
check=[y,y2']              %y==y2但y为行向量，y2为列向量

高阶微分方程 $y^{(n)}=f(t,y,y',\cdots,y^{(n-1)})$ 和一阶微分方程组的解法
高阶微分方程，必须做变量替换，化为一阶微分方程组，才能用MATLAB求数值解；一阶微分方程组的解法和简单的一阶方程解法是一样的。
例6.8：
$y^{'''} - 3 y^{''} - y^{'} y = 0, y (0) = 0, y^{'} (0) = 1, y^{''} (0) = - 1$
解：（1）令 $y 1 = y, y 2 = y^{'}, y 3 = y^{''}$ ，则：
$\left\{ \begin{aligned} y2=y1' \quad \text{y1(0)=0} \\ y2'=y3 \quad \text{y2(0)=1} \\ y3'=3y3+y1y2\quad \text{y3(0)= -1} \end{aligned} \right.$
（2）把一阶方程组写成接受两个参数 $t$ 和 $y$ ,尽管不一定用到参数 $t$ ,必须写两参数。
例6.9

dy=@(t,y)[y(2);y(3);3*y(3)+y(2)*y(1)];
[t,y]=ode15s(dy,[1,3],[0;1;1]);
% y与时刻t是一一对应的，
%y(:,1)是初值问题的解，
%y(:,2)是解的导数，
%y(:,3)是解的二阶导数

plot(t,y(:,1),'*')
title('Solution of van der pol education');   
xlabel('timt t');
ylabel('solution y')

例6.10

rho=10;
beta=28;
lamda=8/3;
f=@(t,y)[rho*(y(2)-y(1));beta*y(1)-y(2)-y(1)*y(3);y(1)*y(2)-lamda*y(3)];
[t,y]=ode45(f,[0,10],[5;13;17])
subplot(2,2,1)
plot(t,y(:,1),'*');
subplot(2,2,2)
%subplot(m,n,p)：m表示是图排成m行，n表示图排成n列，p是指你现在要把曲线画到figure中哪个图上

plot(t,y(:,2),'*');
subplot(2,2,3)
plot(t,y(:,3),'*');

边值问题的MATLAB数值解

MATLAB关于bvp4c命令的官方帮助文档（中文版）：(https://ww2.mathworks.cn/help/matlab/ref/bvp4c.html)
MAATLAB中用bvp4c（需提供猜测的初始值）和bvpinit命令求解常微分方程的两点边值问题，微分方程的标准形式为：
$y'=f(x,y),\quad bc(y(a),y(b))==0$
或：
$\quad bc(y(a),y(b),p)==0$
式中： $p$ 为有关参数； $y, f$ 为向量函数，求解区间为 $[a, b]$ ; $b c$ 为边值条件
$b v p 4 c$ 的调用格式：
$sol=bvp4c(@odeefun,@bcfun,solinit,options,p1,p2,\cdots)$
其中：
（1）函数 $o d e f u n$ 的格式：
$yprime=odefun(x,y,p1,p2,\cdots);$
（2）函数 $b c f u n$ 的格式：
$res=bcfun(ya,yb,p1,p2,\cdots);$
（3）输出参数 $s o l$ 是包含数值解的一个结构，其中 $s o l . x$ 给出了计算数值解的 $x$ 点， $s o l . x (i)$ 处的数值解由 $s o l . y (:, i)$ 给出，类似的， $s o l . x$ 处数值解的一阶导数值由 $s o l . y p (:, i)$ 给出.

例6.11

yprime=@(x,y) [y(2);(y(1)-1)*(1+y(2)^2)^(3/2)];
    %定义一阶方程组的匿名函数

res=@(ya,yb)[ya(1);yb(1)];  
    % ya 和 yb 是与 bc(y(a),y(b),p)=0中的y(a) 和 y(b) 
    %定义边值条件的匿名函数

yinit=@(x)[x.^2;2*x]; 
 %定义初始猜测解的匿名函数，使用x^2与2x作为初始猜测解（初始解是可以任意取的）

solinit=bvpinit(linspace(-1,1,20),yinit);
%linspace(x1,x2,N)，用于产生x1,x2之间的N点行线性的矢量。其中x1、x2、N分别为起始值、终止值、元素个数。若默认N，默认点数为100。
%solinit = bvpinit(x,yinit) 得出边界值问题求解器的初始估计值。x 是指定初始网格的向量。如果要在 [a,b] 中对 BVP 求解，请将 x(1) 指定为 a，并将 x(end) 指定为 b

sol=bvp4c(yprime,res,solinit);

fill(sol.x,sol.y(1,:),[0.7,0.7,0.7])
 %填充解曲线
axis([-1 1 0 1]) 
%axis 设置坐标轴范围和纵横比,axis[x1 x2 x3 x4]表明x轴距离为从x1到x2,y轴距离为从x3到x4

xlabel('x','FontSize',12)  
%'FontSize',12为将字体大小设置为 12 磅

ylabel('h','Rotation',0,'FontSize',12)

例6.12

yprime=@(x,y,mu)[y(2);-mu*y(1)];   %mu为未知参数
res=@(ya,yb,mu)[ya(1);ya(2)-1;yb(1)+yb(2)];  % ya 和 yb 是与 bc(y(a),y(b),p)=0中的y(a) 和 y(b) 
yinit=@(x)[sin(x);cos(x)];
solint=bvpinit(linspace(0,1,10),yinit,5);
sol=bvp4c(yprime,res,solint);
plot(sol.x,sol.y(1,:),'-',sol.x,sol.yp(1,:),'--','LineWidth',2)
xlabel('x','FontSize',12)
legend('y_1','y_2')

例6.13

yprime=@(x,y)[0.5*y(1)*(y(3)-y(1))/y(2);-0.5*(y(3)-y(1))/y(2);(0.9-1000*(y(3)-y(5))-0.5*y(3)*(y(3)-y(1)))/y(4);0.5*(y(3)-y(1));-100*(y(5)-y(3))];
res=@(ya,yb)[ya(1)-1;ya(2)-1;ya(3)-1;ya(4)+10;yb(3)-yb(5)];
yinit=@(x)[1;1;-4.5*x^2+8.91*x+1;-10;-4.5*x^2+9*x+0.91];
solinit=bvpinit(linspace(0,1,5),yinit);
sol=bvp4c(yprime,res,solinit);
plot(sol.x,sol.y(1,:),'-*',sol.x,sol.y(2,:),'-D',sol.x,sol.y(3,:),':s',sol.x,sol.y(4,:),'-.O',sol.x,sol.y(5,:),'--p')
legend('u','v','w','z','y','Location','southwest')    %图注标注在左下角