数学建模暑期集训5：matlab求解常微分方程/偏微分方程

本篇将介绍用matlab求解常微分方程的数值解和解析解，并非是一种完整的模型，仅仅是一些算法。由于数学原理过于复杂，故不探究背后的数学原理，仅将matlab求解的相关函数加以记录。所有代码均可跑通。

1.Matlab求常微分方程的数值解

1.1非刚性常微分方程的数值解法：

功能函数：ode45，ode23，ode113
例：用RK方法（四阶龙格—库塔方法）求解方程
f=-2y+2x^2+2*x

matlab程序：

//doty.m
function f=doty(x,y)
f=-2*y+2*x^2+2*x;
end

//main.m
[x,y]=ode45('doty',[0,0.5],1)

注：[0,0.5]表示求解区间；1为初值列向量

1.2刚性常微分方程的数值解法

功能函数：如ode15s，ode23s，ode23t， ode23tb
使用方法与非刚性类似

1.3高阶微分方程的解法

2.Matlab求常微分方程的解析解

2.1求常微分方程的通解

syms x y
diff_equ='x^2+y+(x-2*y)*Dy=0'
dsolve(diff_equ,'x')

注：'x’代表x为自变量，D代表求导

2.2求常微分方程的初边值问题

syms x y
diff_equ='D3y-D2y=x'
dsolve(diff_equ,'y(1)=8,Dy(1)=7,D2y(2)=4','x')

2.3求常微分方程组

equ1='D2f+3*g=sin(x)';
equ2='Dg+Df=cos(x)';
[general_f,general_g]=dsolve(equ1,equ2,'x')
[f,g]=dsolve(equ1,equ2,'Df(2)=0,f(3)=3,g(5)=1','x')

3.Matlab求解偏微分方程

%(1)问题定义
g='circleg'; %单位圆
b='circleb1'; %边界上为零条件
c=1;a=0;f=1; %（2）产生初始的三角形网格 
[p,e,t]=initmesh(g); %（3）迭代直至得到误差允许范围内的合格解
error=[]; err=1;
while err > 0.01, 
    [p,e,t]=refinemesh(g,p,e,t);
    u=assempde(b,p,e,t,c,a,f); %求得数值解
    exact=(1-p(1,:).^2-p(2,:).^2)/4;
    err=norm(u-exact',inf);
    error=[error err];
end %结果显示
subplot(2,2,1),pdemesh(p,e,t);
subplot(2,2,2),pdesurf(p,t,u)
subplot(2,2,3),pdesurf(p,t,u-exact')

4.Matlab pdetool工具箱求解偏微分方程

对于一般的区域，任意边界条件的偏微分方程，我们可以利用Matlab中pdetool提供的偏微分方程用户图形界面解法。pdetool提供的用户图形界面解法的使用步骤如下：
（i）在Matlab命令窗口运行pdetool，出现PDE Toolbox界面。
（ii）用鼠标点一下工具栏上的“PDE"按钮，在弹出的对话框中定义偏微分方程。
（iii）用鼠标点一下工具栏上的区域按钮，在下面的坐标系中画出偏微分方程的大致定解区域。
（iv）双击（iii）中画出的大致区域，在弹出的对话框中精确定位定解区域。
（v）用鼠标点一下工具栏上的边界按钮“ ”，画出区域的边界。
（vi）双击坐标系中的区域边界，定义偏微分方程的边界条件。
（vii）用鼠标点工具栏上的剖分按钮，对求解区域进行剖分。
（viii）如果求抛物型或双曲型方程的数值解，还需要通过“solve”菜单下的“parameters…”选项设置初值条件。
（ix）用鼠标点一下工具栏上的“＝”按钮，就画出偏微分方程数值解的图形。通过“solve”菜单下的“Export Solution…”选项可以把数值解u输出到Matlab的工作间。
（x）如要画出数值解的三维图形，需要设置“plot"菜单下的“parameters…”选项。

详细操作见
Matlab偏微分方程快速上手：使用pde有限元工具箱求解二维偏微分方程
偏微分方程的数值解(六)： 偏微分方程的 pdetool 解法