偏微分方程的matlab解法
第三步：选取边界
首先选择Boundary菜单中Boundary Mode 命令，进入边界模式．然后单击Boundary菜单 中Remove All Subdomain Borders选项。从而 去掉子域边界，如图22 2．单击Boundary菜单 中Specify Boundary Conditions选项，打开 Boundary Conditions对话框，输入边界件．本 例取默认条件，即将全部边界设为齐次Dirichlet 条件，边界显示为红色． 如果想将几何与边界信息存储，可选Boundary 菜单中的Export Decomposed Geometrv．Boundary Cond's命令，将它们分 别存储在g、b变量中，并通过MATLAB形成M文 件．
求解双曲型方程的例子
例24.2.1 用 MATLAB 求解下面波动方程定解问题并动态显示解的分布
2u (2u t 2 x2
2u ) 0 y 2
u
|x
1
u
|x1
0,
u y
y 1u Fra bibliotekyy1 0
π
π
u(x,
y, 0)
atan[ sin(
保持在100 °C，板的右边热量从板向环境空气定常流动，
t t 其他边及内孔边界保持绝缘。初始
°C ，于是概括为如下定解问题;
是板的温度为0 0
d u u0, t
u 100 ,在左边界上
u 1，在右边界上 n u = 0，其他边界上 n
u tto 0
区域的边界顶点坐标为(-0.5,-0.8), (0.5,-0.8), (-0.5,0.8), (0.5,0.8)。 内边界顶点坐标(-0.05,-0.4), (-0.05,0.4) ,(0.05,-0.4), (0.05,0.4)。
第五步：选择Mesh菜单中Initialize Mesh命令，
进行网格剖分， 选择Mesh菜单中Refine Mesh命令，使网格密集化，
如图22.3．
图 22.3 网格密集化
第六步： 解偏微分方程并显示图形解
选择Solve菜单中Solve PDE命令，解 偏微分方程并显示图形解，如图 2.4 所示
例： 解热传导方程 ut u f
边界条件是齐次类型，定解区域自定。
【解】 第一步：启动MATLAB，键入命令pdetool并回车，就 进入GUI．在Options菜单下选择Gid命令，打开栅格， 栅格使用户容易确定所绘图形的大小． 第二步：选定定解区域本题为自定区域：自拟定解区 域如图22 1所示：E1-E2+R1-E3．具体用快捷工具分别 画椭圆E1、圆E2、矩形R1、圆E3．然后在Set formula 栏中进行编辑并用算术运算符将图形对象名称连接起 来(或删去默认的表达式，直接键入E1-E2+R1-E3)
例如，对于细杆导热，虽然是一维问题， 可以将宽度y虚拟出来，对应于y的边界 条件和初始条件按照题意制定
Boundary Mode
PDE Mode
PDE Specification，确定偏 微分方程类型共有四种：
椭圆形Elliptic
抛物型Parabolic
双曲型Hyperbolic
第七步：单击Plot菜单中Parameter选项，打开Plot Selection对话框，选中Color,Height(3D plot)和 Show mesh三项.再单击Polt按钮，显示三维图形解， 如图22.5所示.
第八步：若要画等值线图和矢量场图，单击plot菜单 中parameter 选项，在plot selection对话框中选中 contour 和arrow两选项。然后单击plot按钮，可显 示解的等值线图和矢量场图，如图2.6所示。
图 2.6 解的等值线图和矢量场图
求解椭圆型方程的例子
单位圆上的poisson方程边值问题：
-u=1 , = (x, y) x2 y2 1 ,
u 0 问题的精确解为
u(x,y)= (1 x2 y2 ) . 4
求解抛物型方程的例子
考虑一个带有矩形孔的金属板上的热传导问题。板的左边
3、解的可视化．
PDEToolbox注意事项
只能解决二维模型，一维的扩成二维，三 维的缩成二维，时间维不计算在内 公式类型，只能解决部分偏微分方程，由 公式类型决定 边界条件两种，Dirichlet和Neumann 初始条件
先确定方程大类
Draw Mode
画图模式，先将处理的区域画出来，二 维，方形，圆形，支持多边形，可以手 动更改坐标，旋转rotate
图 22.2 定解问题的边界
第四步：设置方程类型
选择PDE菜单中PDE Mode命令，进入PDE模式， 再单击PDE菜单中PDE Secification选项，打开 PDE Secification对话框，设置方程类型．
本例取抛物型方程 d u (cu) au f ,
t
故参数c,a,f,d,分别是l，0，10，1．
2
x)], ut
( x,
y,
0)
2
cos(πx)
exp[cos(
2
y)]
已知求解域是方形区域，空间坐标的个数由具体问题确定．
偏微分方程的matlab解法
主要讲述如何用MATLAB实现对偏微分方程的仿 真求解．MATLAB的偏微分方程工具箱（PDE Toolbox）的出现，为偏微分方程的求解以及定 性研究提供了捷径．主要步骤为：
1、设置PDE的定解问题．即设置二维定解区域、 边界条件以及方程的形式；和系数
2、用有限元法（FEM）求解PDE．即网格的生 成、方程的离散以及求出数值解；
Mesh Mode
网格划分，细化
Solve，Plot
如果有初始条件（与t有关），则在 Solve的Parameters里有其设定，如果 没有初始条件（与t无关），则不必设 定Plot只是确定画图的参数，包括是否 动画，是否3D，是否画出等温线，是否 有箭头。。。
Save As
保存成M-file，自动生成
偏微分方程Matlab数值解法(补充4)
偏微分方程Matlab 数值解法（补充4）Matlab 可以求解一般的偏微分方程，也可以利用偏微分方程工具箱中给出的函数求解一些偏微分方程。
1 偏微分方程组求解Matlab 语言提供了pdepe()函数，可以直接求解偏微分方程(,,,)[(,,,)](,,,)m mu u u u C x t u x x f x t u s x t u x t x x x-∂∂∂∂∂=+∂∂∂∂∂ （4.1）这样，偏微分方程可以编写为以下函数的描述，其入口为[,,](,,,)x c f s pdefun x t u u =其中：pdefun 为函数名。
由给定输入变量可计算出,,c f s 这三个函数。
边界条件可以用下面的函数描述(,,)(,,).*(,,,)0up x t u q x t u f x t u x∂+=∂ （4.2） 这样的边值函数可以写为Matlab 函数[,,,](,,,)a a b b x p q p q pdebc x t u u =初始条件数学描述为00(,)u x t u =，编写一个简单的函数即可0()u pdeic x =还可以选择x 和t 的向量，再加上描述这些函数，就可以用pdepe （）函数求解次偏微分方程，需要用如下格式求解(,@,@,@,,)sol pdepe m pdefun pdeic pdebc x t =【例1】 试求下列偏微分方程2111222221220.024()0.17()u u F u u t x u u F u u tx ⎧∂∂=--⎪⎪∂∂⎨∂∂⎪=+-⎪∂∂⎩ 其中： 5.7311.46()xx F x e e -=-，且满足初始条件1(,0)1u x =，2(,1)0u x =及边界条件1221(0,)0,(0,)0,(1,)1,(1,)0u u t u t u t t x x∂∂====∂∂解：对照给出的偏微分方程和（4.1），可将原方程改写为111222120.024/()1.*10.17/()u u x F u u u u x F u u t x ∂∂--⎡⎤⎡⎤⎡⎤⎡⎤∂∂=+⎢⎥⎢⎥⎢⎥⎢⎥∂∂-∂∂⎣⎦⎣⎦⎣⎦⎣⎦可知0m =，且1122120.024/()1,,10.17/()u x F u u c f s u x F u u ∂∂--⎡⎤⎡⎤⎡⎤===⎢⎥⎢⎥⎢⎥∂∂-⎣⎦⎣⎦⎣⎦编写下面的Matlab 函数function [c,f,s]=c7mpde(x,t,u,du)c=[1;1];y=u(1)-u(2);F=exp(5.73*y)-exp(-11.46*y);s=F*[-1;1]; f=[0.024*du(1);0.17*du(2)];套用（4.2）中的边界条件，可以写出如下的边值方程左边界2010.*00f u ⎡⎤⎡⎤⎡⎤+=⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦，右边界1100.*100u f -⎡⎤⎡⎤⎡⎤+=⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦ 编写下面的Matlab 函数function [pa,qa,pb,qb]=c7mpbc(xa,ua,xb,ub,t) pa=[0;ua(2)];qa=[1;0];pb=[ub(1)-1;0];qb=[0,1]; 另外，描述初值的函数function u0=c7mpic(x) u0=[1;0];有了这三个函数，选定x 和t 向量，则可以由下面的程序直接求此微分方程，得出解1u 和2u 。
matlab求解最简单的一阶偏微分方程
matlab求解最简单的一阶偏微分方程一、引言在科学和工程领域，偏微分方程是非常重要的数学工具，用于描述各种现象和过程。
而MATLAB作为一种强大的数值计算软件，可以用来求解各种复杂的偏微分方程。
本文将以MATLAB求解最简单的一阶偏微分方程为主题，探讨其基本原理、数值求解方法以及具体实现过程。
二、一阶偏微分方程的基本原理一阶偏微分方程是指只含有一个未知函数的偏导数的微分方程。
最简单的一阶偏微分方程可以写成如下形式：\[ \frac{\partial u}{\partial t} = F(x, t, u, \frac{\partial u}{\partial x}) \]其中，\(u(x, t)\) 是未知函数，\(F(x, t, u, \frac{\partial u}{\partial x})\) 是给定的函数。
一阶偏微分方程可以描述很多实际问题，比如热传导、扩散等。
在MATLAB中，我们可以使用数值方法求解这类方程。
三、数值求解方法1. 有限差分法有限差分法是一种常用的数值求解偏微分方程的方法。
其基本思想是用离散的方式来逼近偏导数，然后将偏微分方程转化为代数方程组。
在MATLAB中，我们可以使用内置的求解器来求解离散化后的代数方程组。
2. 特征线法特征线法是另一种常用的数值求解方法，它利用特征线方程的特点来求解偏微分方程。
这种方法在求解一维情况下的偏微分方程时特别有效，可以提高求解的效率和精度。
四、MATLAB求解过程在MATLAB中，我们可以使用`pdepe`函数来求解一阶偏微分方程。
该函数可以针对特定的方程和边界条件，利用有限差分法进行离散化求解。
下面给出一个具体的例子来说明如何使用MATLAB求解最简单的一阶偏微分方程。
假设我们要求解如下的一维热传导方程：\[ \frac{\partial u}{\partial t} = \alpha \frac{\partial^2 u}{\partial x^2} \]其中，\(\alpha\) 是热传导系数。
(完整版)偏微分方程的MATLAB解法
引言偏微分方程定解问题有着广泛的应用背景。
人们用偏微分方程来描述、解释或者预见各种自然现象，并用于科学和工程技术的各个领域fll。
然而，对于广大应用工作者来说，从偏微分方程模型出发，使用有限元法或有限差分法求解都要耗费很大的工作量，才能得到数值解。
现在，MATLAB PDEToolbox已实现对于空间二维问题高速、准确的求解过程。
偏微分方程如果一个微分方程中出现的未知函数只含一个自变量，这个方程叫做常微分方程，也简称微分方程；如果一个微分方程中出现多元函数的偏导数，或者说如果未知函数和几个变量有关，而且方程中出现未知函数对几个变量的导数，那么这种微分方程就是偏微分方程。
常用的方法有变分法和有限差分法。
变分法是把定解问题转化成变分问题，再求变分问题的近似解；有限差分法是把定解问题转化成代数方程，然后用计算机进行计算；还有一种更有意义的模拟法，它用另一个物理的问题实验研究来代替所研究某个物理问题的定解。
虽然物理现象本质不同，但是抽象地表示在数学上是同一个定解问题，如研究某个不规则形状的物体里的稳定温度分布问题，由于求解比较困难，可作相应的静电场或稳恒电流场实验研究，测定场中各处的电势，从而也解决了所研究的稳定温度场中的温度分布问题。
随着物理科学所研究的现象在广度和深度两方面的扩展，偏微分方程的应用范围更广泛。
从数学自身的角度看，偏微分方程的求解促使数学在函数论、变分法、级数展开、常微分方程、代数、微分几何等各方面进行发展。
从这个角度说，偏微分方程变成了数学的中心。
一、MATLAB方法简介及应用1.1 MATLAB简介MATLAB是美国MathWorks公司出品的商业数学软件，用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境，主要包括MATLAB和Simulink两大部分。
1.2 Matlab主要功能数值分析数值和符号计算工程与科学绘图控制系统的设计与仿真数字图像处理数字信号处理通讯系统设计与仿真财务与金融工程1.3 优势特点1) 高效的数值计算及符号计算功能，能使用户从繁杂的数学运算分析中解脱出来；2) 具有完备的图形处理功能，实现计算结果和编程的可视化；3) 友好的用户界面及接近数学表达式的自然化语言，使学者易于学习和掌握；4) 功能丰富的应用工具箱(如信号处理工具箱、通信工具箱等) ，为用户提供了大量方便实用的处理工具。
matlab解二阶偏微分方程
<div align="center">一、Matlab求解二阶偏微分方程(ODE)的基本步骤</div>1. 数学模型：首先要确定求解的方程是哪一类的偏微分方程(ODE)，然后建立其对应的数学模型，使其符合这类微分方程的形式;2. 确定边界条件：确定迭代范围$[a,b]$,边界条件函数 $y(a)=\alpha$ 、$y(b)=\beta$;3. 写出Matlab程序：在该类ODE中，通常会有某一种常用的数值求解方法，一般使用微分方程求解器(ODE)，如ode45等;4. 获得实际结果：开始编写Matlab程序，完成参差和参数的输入以后，可以运行Matlab程序，然后求得结果，再用图像表示出来。
<div align="center">二、具体求解</div>$$\frac{d^2y}{dx^2}+y=6sin(2x)$$微分方程为二阶常微分方程，求解条件如下：$[a,b]=[0,\pi], y(0)=1,y(\pi)=3.$（1）Matlab函数表达式首先建立与二阶非齐次线性常微分方程相符合的数学模型，其Matlab函数表达式为$$ f(x,y,y')=\frac{dy}{dx}-y'-6sin2x $$其中，$y=y(x)$；（2）函数程序在Matlab中，定义函数程序 $myode.m$ ，此程序返回右端函数 $f(x,y,y')$ 的值表达式，程序内容如下。
```MATLAB% 右端函数程序function dy=myode(x,y)dy=[y(2);-y(2)-6*sin(2*x)];end```（3）调用Matlab函数olvede45调用Matlab函数 solvede45 求解二阶ODN，程序内容如下：```MATLAB% 主程序求解% maxstep表示分裂的步长大小% Tolerence表示误差，控制求解精度Maxstep=0.25;Tolerence=1e-4;a=0;b=pi;y0=[1;0];[x,y] = ode45('myode',[a,b],y0,options);```（4）结果展示输入参数之后，运行Matlab程序，得到如下图：此图为$y(x)$随$x$变化的曲线，可以看出，二阶偏微分的求解结果满足了边界条件，即$y(0)=1,y(\pi)=3$ ，如图中红色圆点所示。
(完整版)偏微分方程的MATLAB解法
引言偏微分方程定解问题有着广泛的应用背景。
人们用偏微分方程来描述、解释或者预见各种自然现象，并用于科学和工程技术的各个领域fll。
然而，对于广大应用工作者来说，从偏微分方程模型出发，使用有限元法或有限差分法求解都要耗费很大的工作量，才能得到数值解。
现在，MATLAB PDEToolbox已实现对于空间二维问题高速、准确的求解过程。
偏微分方程如果一个微分方程中出现的未知函数只含一个自变量，这个方程叫做常微分方程，也简称微分方程；如果一个微分方程中出现多元函数的偏导数，或者说如果未知函数和几个变量有关，而且方程中出现未知函数对几个变量的导数，那么这种微分方程就是偏微分方程。
常用的方法有变分法和有限差分法。
变分法是把定解问题转化成变分问题，再求变分问题的近似解；有限差分法是把定解问题转化成代数方程，然后用计算机进行计算；还有一种更有意义的模拟法，它用另一个物理的问题实验研究来代替所研究某个物理问题的定解。
虽然物理现象本质不同，但是抽象地表示在数学上是同一个定解问题，如研究某个不规则形状的物体里的稳定温度分布问题，由于求解比较困难，可作相应的静电场或稳恒电流场实验研究，测定场中各处的电势，从而也解决了所研究的稳定温度场中的温度分布问题。
随着物理科学所研究的现象在广度和深度两方面的扩展，偏微分方程的应用范围更广泛。
从数学自身的角度看，偏微分方程的求解促使数学在函数论、变分法、级数展开、常微分方程、代数、微分几何等各方面进行发展。
从这个角度说，偏微分方程变成了数学的中心。
一、MATLAB方法简介及应用1.1 MATLAB简介MATLAB是美国MathWorks公司出品的商业数学软件，用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境，主要包括MATLAB和Simulink两大部分。
1.2 Matlab主要功能数值分析数值和符号计算工程与科学绘图控制系统的设计与仿真数字图像处理数字信号处理通讯系统设计与仿真财务与金融工程1.3 优势特点1) 高效的数值计算及符号计算功能，能使用户从繁杂的数学运算分析中解脱出来；2) 具有完备的图形处理功能，实现计算结果和编程的可视化；3) 友好的用户界面及接近数学表达式的自然化语言，使学者易于学习和掌握；4) 功能丰富的应用工具箱(如信号处理工具箱、通信工具箱等) ，为用户提供了大量方便实用的处理工具。
matlab解偏微分方程
ui,j +1 − ui,j = Hui,j +1 ∆t Hui,j = a2 ui+1,j − 2ui,j + ui−1,j (∆x)2
ui,j +1 − ui,j = Hui,j ∆t 将显式与隐式相加,得平均公式 ui,j +1 − ui,j 1 1 = Hui,j + Hui,j +1 ∆t 2 2
得ui,0 = ui,2 − 2ψi
1 ui,2 = [c(ui+1,1 + ui−1,1) + 2(1 − c)ui,1 + 2ψi t] 2
3.3
例题 两端固定的弦振动
两端固定的弦, 初速为零，初位移是 h x, (0 ≤ x ≤ 2/3) 2 / 3 u(x, 0) = 1−x , (2/3 < x ≤ 1) h 1 − 2/3
作图所用程序如下,其中取c = 0.05, l = 1, h = 0.05.这里使用的方程 与初始条件表示方法与上一节相同. N=4000; c=0.05; x=linspace(0,1,420)’; u1(1:420)=0; u2(1:420)=0; u3(1:420)=0; u1(2:280)=0.05/279*(1:279)’; u1(281:419)=0.05/(419-281)*(419-(281:419)’); u2(2:419)=u1(2:419)+c/2*(u1(3:420)-2*u1(2:419)+u1(1:418)); h=plot(x,u1,’linewidth’,3); axis([0,1,-0.05,0.05]); set(h,’EraseMode’,’xor’,’MarkerSize’,18) for k=2:N set(h,’XData’,x,’YData’,u2) ; drawnow; u3(2:419)=2*u2(2:419)-u1(2:419)+c*(u2(3:420)... -2*u2(2:419)+u2(1:418)); u1=u2; u2=u3; end
【精品】偏微分的MATLAB数值解法课件
方法一：pdepe函数实现
• x=0:1:40; • t=0:0.01:0.2; • m=0; • sol=pdepe(m,@pdefun,@pdeic,@pdebc,x,t); • b=sol(20,:); • plot(x,b); • title('the solution of u') • xlabel('x') • ylabel('y') • zlabel('u')
偏微分的MATLAB数值解法
偏微分的MATLAB数值解法
• 方法一：pdepe函数实现 • 方法二：pdetool实现 • 方法三：程序实现
方法一：pdepe函数实现
• @pdeic： • function u0=pdeic(x) • if x<10 • u0=0; • elseif x<30 • u0=1; • else • u0=0; • end
•
end
• end
方法三：程序实现
图 22.12 波动方程解析解的分布
偏微分的MATLAB数值解法
• 方法总结： • 1.pdede调用简单，但计算功能稍弱 • 2.pdetool使用方便，但限于四种方程类
型 • 3.程序编写较为繁琐
方法一：pdepe函数实现
方法二：pdetool实现
• 1.pdetool界面 • 2.选定求解微分方程类型（双曲线、抛物线、椭
圆、特殊值型）并设定参数 • 3.绘制求解区域 • 4.边界条件和初值条件（Dirichlet和Neumann） • 5.生成网格 • 6.求解方程并绘制图形
方法二：pdetool实现
• 应用实例：
u(ux,y)
x2 y x0
matlab 偏微分方程
MATLAB是一个强大的数值计算环境，可以用来解决各种各样的数学问题，包括偏微分方程。
下面是一个简单的例子，展示如何在MATLAB中解决一维的偏微分方程。
假设我们要解决以下一维的热传导方程：
∂u∂t=∂2u∂x2
在给定的初始条件和边界条件下：
u(x,0)=sin(πx)u(0,t)=0, u(1,t)=0
我们可以使用MATLAB中的pdepe函数来求解这个问题。
以下是一个简单的MATLAB代码示例：
```matlab
定义参数
T = 1; 最终时间
h = 0.01; 空间步长
t = 0:T/h:T; 时间向量
x = 0:h:1; 空间向量
n = length(x); 空间点的数量
m = length(t); 时间点的数量
初始化矩阵存储解
U = zeros(m, n);
U(:,1) = sin(pi*x); 初始条件
定义偏微分方程
pdepe('u_tt', U, t, x, 'heat', 'periodic');
使用pdepe求解偏微分方程
[U, ~] = pdepe(U, t, x);
绘制结果
surf(x, t, U);
```
这个代码示例使用了MATLAB的pdepe函数，这是一个用于求解偏微分方程的函数。
在上面的代码中，我们首先定义了参数，然后初始化了存储解的矩阵。
然后，我们定义了偏微分方程，并使用pdepe 函数求解它。
最后，我们使用surf函数绘制了结果。
matlab求解初边值问题的偏微分方程
偏微分方程是描述自然界中动态过程的重要数学工具，在工程领域中，求解偏微分方程是很多实际问题的重要一步。
MATLAB作为一种强大的数学计算软件，提供了丰富的工具和函数来求解各种类型的偏微分方程。
本文将介绍使用MATLAB求解初边值问题的偏微分方程的方法和步骤。
一、MATLAB中的偏微分方程求解工具MATLAB提供了几种可以用来求解偏微分方程的工具和函数，主要包括：1. pdepe函数：用于求解偏微分方程初边值问题的函数，可以处理各种类型的偏微分方程，并且可以自定义边界条件和初始条件。
2. pdepeplot函数：用于绘制pdepe函数求解得到的偏微分方程的解的可视化图形，有助于直观地理解方程的解的特性。
3. pdetool工具箱：提供了一个交互式的图形用户界面，可以用来建立偏微分方程模型并进行求解，适用于一些复杂的偏微分方程求解问题。
二、使用pdepe函数求解偏微分方程初边值问题的步骤对于给定的偏微分方程初边值问题，可以按照以下步骤使用pdepe函数进行求解：1. 定义偏微分方程需要将给定的偏微分方程转化为标准形式，即将偏微分方程化为形式为c(x,t,u)∂u/∂t = x(r ∂u/∂x) + ∂(p∂u/∂x) + f(x,t,u)的形式。
2. 编写边界条件和初始条件函数根据实际问题的边界条件和初始条件，编写相应的函数来描述这些条件。
3. 设置空间网格选择合适的空间网格来离散空间变量，可以使用linspace函数来产生均匀分布的网格。
4. 调用pdepe函数求解偏微分方程将定义好的偏微分方程、边界条件和初始条件函数以及空间网格作为参数传递给pdepe函数，调用该函数求解偏微分方程。
5. 可视化结果使用pdepeplot函数绘制偏微分方程的解的可视化图形，以便对解的性质进行分析和理解。
三、实例分析考虑一维热传导方程初边值问题：∂u/∂t = ∂^2u/∂x^2, 0<x<1, 0<t<1u(0,t) = 0, u(1,t) = 0, u(x,0) = sin(πx)使用MATLAB求解该初边值问题的步骤如下：1. 定义偏微分方程将热传导方程化为标准形式，得到c(x,t,u) = 1, r = 1, p = 1, f(x,t,u) = 0。
Matlab偏微分方程求解方法
Matlab 偏微分方程求解方法目录：§1 Function Summary on page 10-87§2 Initial Value Problems on page 10-88§3 PDE Solver on page 10-89§4 Integrator Options on page 10-92§5 Examples” on page 10-93§1 Function Summary1.1 PDE Solver” on page 10-871,2 PDE Helper Functi on” on page 10-871.3 PDE SolverThis is the MATLAB PDE solver.PDE Helper Function§2 Initial Value Problemspdepe solves systems of parabolic and elliptic PDEs in one spatial variable x and time t, of the form)xu ,u ,t ,x (s ))x u ,u ,t ,x (f x (x x t u )x u ,u ,t ,x (c m m ∂∂+∂∂∂∂=∂∂∂∂- (10-2) The PDEs hold for b x a ,t t t f 0≤≤≤≤.The interval [a, b] must be finite. mcan be 0, 1, or 2, corresponding to slab, cylindrical, or spherical symmetry,respectively. If m > 0, thena ≥0 must also hold.In Equation 10-2,)x /u ,u ,t ,x (f ∂∂ is a flux term and )x /u ,u ,t ,x (s ∂∂ is a source term. The flux term must depend on x /u ∂∂. The coupling of the partial derivatives with respect to time is restricted to multiplication by a diagonal matrix )x /u ,u ,t ,x (c ∂∂. The diagonal elements of this matrix are either identically zero or positive. An element that is identically zero corresponds to an elliptic equation and otherwise to a parabolic equation. There must be at least one parabolic equation. An element of c that corresponds to a parabolic equation can vanish at isolated values of x if they are mesh points.Discontinuities in c and/or s due to material interfaces are permitted provided that a mesh point is placed at each interface.At the initial time t = t0, for all x the solution components satisfy initial conditions of the form)x (u )t ,x (u 00= (10-3)At the boundary x = a or x = b, for all t the solution components satisfy a boundary condition of the form0)xu ,u ,t ,x (f )t ,x (q )u ,t ,x (p =∂∂+ (10-4) q(x, t) is a diagonal matrix with elements that are either identically zero or never zero. Note that the boundary conditions are expressed in terms of the f rather than partial derivative of u with respect to x-x /u ∂∂. Also, ofthe two coefficients, only p can depend on u.§3 PDE Solver3.1 The PDE SolverThe MATLAB PDE solver, pdepe, solves initial-boundary value problems for systems of parabolic and elliptic PDEs in the one space variable x and time t.There must be at least one parabolic equation in the system.The pdepe solver converts the PDEs to ODEs using a second-order accurate spatial discretization based on a fixed set of user-specified nodes. The discretization method is described in [9]. The time integration is done with ode15s. The pdepe solver exploits the capabilities of ode15s for solving the differential-algebraic equations that arise when Equation 10-2 contains elliptic equations, and for handling Jacobians with a specified sparsity pattern. ode15s changes both the time step and the formula dynamically.After discretization, elliptic equations give rise to algebraic equations. If the elements of the initial conditions vector that correspond to elliptic equations are not “consistent” with the discretization, pdepe tries to adjust them before eginning the time integration. For this reason, the solution returned for the initial time may have a discretization error comparable to that at any other time. If the mesh is sufficiently fine, pdepe can find consistent initial conditions close to the given ones. If pdepe displays amessage that it has difficulty finding consistent initial conditions, try refining the mesh. No adjustment is necessary for elements of the initial conditions vector that correspond to parabolic equations.PDE Solver SyntaxThe basic syntax of the solver is:sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan)Note Correspondences given are to terms used in “Initial Value Problems” on page 10-88.The input arguments arem: Specifies the symmetry of the problem. m can be 0 =slab, 1 = cylindrical, or 2 = spherical. It corresponds to m in Equation 10-2. pdefun: Function that defines the components of the PDE. Itcomputes the terms f,c and s in Equation 10-2, and has the form[c,f,s] = pdefun(x,t,u,dudx)where x and t are scalars, and u and dudx are vectors that approximate the solution and its partial derivative with respect to . c, f, and s are column vectors. c stores the diagonal elements of the matrix .icfun: Function that evaluates the initial conditions. It has the formu = icfun(x)When called with an argument x, icfun evaluates and returns the initial values of the solution components at x in the column vector u.bcfun:Function that evaluates the terms and of the boundary conditions. Ithas the form[pl,ql,pr,qr] = bcfun(xl,ul,xr,ur,t)where ul is the approximate solution at the left boundary xl = a and ur is the approximate solution at the right boundary xr = b. pl and ql are column vectors corresponding to p and the diagonal of q evaluated at xl. Similarly, pr and qr correspond to xr. When m>0 and a = 0, boundedness of the solution near x = 0 requires that the f vanish at a = 0. pdepe imposes this boundary condition automatically and it ignores values returned in pl and ql.xmesh:Vector [x0, x1, ..., xn] specifying the points at which a numerical solution is requested for every value in tspan. x0 and xn correspond to a and b , respectively. Second-order approximation to the solution is made on the mesh specified in xmesh. Generally, it is best to use closely spaced mesh points where the solution changes rapidly. pdepe does not select the mesh in automatically. You must provide an appropriate fixed mesh in xmesh. The cost depends strongly on the length of xmesh. When , it is not necessary to use a fine mesh near to x=0 account for the coordinate singularity.The elements of xmesh must satisfy x0 < x1 < ... < xn.The length of xmesh must be ≥3.tspan:Vector [t0, t1, ..., tf] specifying the points at which a solution is requested for every value in xmesh. t0 and tf correspond tot and f t,respectively.pdepe performs the time integration with an ODE solver that selects both the time step and formula dynamically. The solutions at the points specified in tspan are obtained using the natural continuous extension of the integration formulas. The elements of tspan merely specify where you want answers and the cost depends weakly on the length of tspan.The elements of tspan must satisfy t0 < t1 < ... < tf.The length of tspan must be ≥3.The output argument sol is a three-dimensional array, such that•sol(:,:,k) approximates component k of the solution .•sol(i,:,k) approximates component k of the solution at time tspan(i) and mesh points xmesh(:).•sol(i,j,k) approximates component k of the solution at time tspan(i) and the mesh point xmesh(j).4.2 PDE Solver OptionsFor more advanced applications, you can also specify as input arguments solver options and additional parameters that are passed to the PDE functions.options:Structure of optional parameters that change the default integration properties. This is the seventh input argument.sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan,options)See “Integrator Options” on page 10-92 for more information.Integrator OptionsThe default integration properties in the MATLAB PDE solver are selected to handle common problems. In some cases, you can improve solver performance by overriding these defaults. You do this by supplying pdepe with one or more property values in an options structure.sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan,options)Use odeset to create the options structure. Only those options of the underlying ODE solver shown in the following table are available for pdepe.The defaults obtained by leaving off the input argument options are generally satisfactory. “Integrator Options” on page 10-9 tells you how to create the structure and describes the properties.PDE Properties§4 Examples•“Single PDE” on page 10-93•“System of PDEs” on page 10-98•“Additional Examples” on page 10-1031.Single PDE• “Solving the Equation” on page 10-93• “Evaluating the Solution” on page 10-98Solving the Equation. This example illustrates the straightforward formulation, solution, and plotting of the solution of a single PDE222x u t u ∂∂=∂∂π This equation holds on an interval 1x 0≤≤ for times t ≥ 0. At 0t = the solution satisfies the initial condition x sin )0,x (u π=.At 0x =and 1x = , the solution satisfies the boundary conditions0)t ,1(xu e ,0)t ,0(u t =∂∂+π=- Note The demo pdex1 contains the complete code for this example. The demo uses subfunctions to place all functions it requires in a single MATLAB file.To run the demo type pdex1 at the command line. See “PDE Solver Syntax” on page 10-89 for more information. 1 Rewrite the PDE. Write the PDE in the form)xu ,u ,t ,x (s ))x u ,u ,t ,x (f x (x x t u )x u ,u ,t ,x (c m m ∂∂+∂∂∂∂=∂∂∂∂- This is the form shown in Equation 10-2 and expected by pdepe. For this example, the resulting equation is0x u x x x t u 002+⎪⎭⎫ ⎝⎛∂∂∂∂=∂∂π with parameter and the terms 0m = and the term0s ,xu f ,c 2=∂∂=π=2 Code the PDE. Once you rewrite the PDE in the form shown above (Equation 10-2) and identify the terms, you can code the PDE in a function that pdepe can use. The function must be of the form[c,f,s] = pdefun(x,t,u,dudx)where c, f, and s correspond to the f ,c and s terms. The code below computes c, f, and s for the example problem.function [c,f,s] = pdex1pde(x,t,u,DuDx)c = pi^2;f = DuDx;s = 0;3 Code the initial conditions function. You must code the initial conditions in a function of the formu = icfun(x)The code below represents the initial conditions in the function pdex1ic. Partial Differential Equationsfunction u0 = pdex1ic(x)u0 = sin(pi*x);4 Code the boundary conditions function. You must also code the boundary conditions in a function of the form[pl,ql,pr,qr] = bcfun(xl,ul,xr,ur,t) The boundary conditions, written in the same form as Equation 10-4, are0x ,0)t ,0(x u .0)t ,0(u ==∂∂+and1x ,0)t ,1(xu .1e t ==∂∂+π- The code below evaluates the components )u ,t ,x (p and )u ,t ,x (q of the boundary conditions in the function pdex1bc.function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t)pl = ul;ql = 0;pr = pi * exp(-t);qr = 1;In the function pdex1bc, pl and ql correspond to the left boundary conditions (x=0 ), and pr and qr correspond to the right boundary condition (x=1).5 Select mesh points for the solution. Before you use the MATLAB PDE solver, you need to specify the mesh points at which you want pdepe to evaluate the solution. Specify the points as vectors t and x.The vectors t and x play different roles in the solver (see “PDE Solver” on page 10-89). In particular, the cost and the accuracy of the solution depend strongly on the length of the vector x. However, the computation is much less sensitive to the values in the vector t.10 CalculusThis example requests the solution on the mesh produced by 20 equally spaced points from the spatial interval [0,1] and five values of t from thetime interval [0,2].x = linspace(0,1,20);t = linspace(0,2,5);6 Apply the PDE solver. The example calls pdepe with m = 0, the functions pdex1pde, pdex1ic, and pdex1bc, and the mesh defined by x and t at which pdepe is to evaluate the solution. The pdepe function returns the numerical solution in a three-dimensional array sol, wheresol(i,j,k) approximates the kth component of the solution,u, evaluated atkt(i) and x(j).m = 0;sol = pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t);This example uses @ to pass pdex1pde, pdex1ic, and pdex1bc as function handles to pdepe.Note See the function_handle (@), func2str, and str2func reference pages, and the @ section of MATLAB Programming Fundamentals for information about function handles.7 View the results. Complete the example by displaying the results:a Extract and display the first solution component. In this example, the solution has only one component, but for illustrative purposes, the example “extracts” it from the three-dimensional array. The surface plot shows the behavior of the solution.u = sol(:,:,1);surf(x,t,u)title('Numerical solution computed with 20 mesh points') xlabel('Distance x') ylabel('Time t')Distance xNumerical solution computed with 20 mesh points.Time tb Display a solution profile at f t , the final value of . In this example,2t t f ==.figure plot(x,u(end,:)) title('Solution at t = 2') xlabel('Distance x') ylabel('u(x,2)')Solutions at t = 2.Distance xu (x ,2)Evaluating the Solution. After obtaining and plotting the solution above, you might be interested in a solution profile for a particular value of t, or the time changes of the solution at a particular point x. The kth column u(:,k)(of the solution extracted in step 7) contains the time history of the solution at x(k). The jth row u(j,:) contains the solution profile at t(j). Using the vectors x and u(j,:), and the helper function pdeval, you can evaluate the solution u and its derivative at any set of points xout [uout,DuoutDx] = pdeval(m,x,u(j,:),xout)The example pdex3 uses pdeval to evaluate the derivative of the solution at xout = 0. See pdeval for details.2. System of PDEsThis example illustrates the solution of a system of partial differential equations. The problem is taken from electrodynamics. It has boundary layers at both ends of the interval, and the solution changes rapidly for small . The PDEs are)u u (F xu017.0t u )u u (F xu 024.0t u 212222212121-+∂∂=∂∂--∂∂=∂∂ where )y 46.11exp()y 73.5exp()y (F --=. The equations hold on an interval1x 0≤≤ for times 0t ≥.The solution satisfies the initial conditions0)0,x (u ,1)0,x (u 21≡≡and boundary conditions0)t ,1(xu,0)t ,1(u ,0)t ,0(u ,0)t ,0(x u 2121=∂∂===∂∂ Note The demo pdex4 contains the complete code for this example. The demo uses subfunctions to place all required functions in a single MATLAB file. To run this example type pdex4 at the command line.1 Rewrite the PDE. In the form expected by pdepe, the equations are⎥⎦⎤⎢⎣⎡---+⎥⎦⎤⎢⎣⎡∂∂∂∂∂∂=⎥⎦⎤⎢⎣⎡∂∂⎥⎦⎤⎢⎣⎡)u u (F )u u (F )x /u 170.0)x /u (024.0x u u t *.1121212121 The boundary conditions on the partial derivatives of have to be written in terms of the flux. In the form expected by pdepe, the left boundary condition is⎥⎦⎤⎢⎣⎡=⎥⎦⎤⎢⎣⎡∂∂∂∂⎥⎦⎤⎢⎣⎡+⎥⎦⎤⎢⎣⎡00)x /u (170.0)x /u (024.0*.01u 0212and the right boundary condition is⎥⎦⎤⎢⎣⎡=⎥⎦⎤⎢⎣⎡∂∂∂∂⎥⎦⎤⎢⎣⎡+⎥⎦⎤⎢⎣⎡-00)x /u (170.0)x /u (024.0*.1001u 2112 Code the PDE. After you rewrite the PDE in the form shown above, you can code it as a function that pdepe can use. The function must be of the form[c,f,s] = pdefun(x,t,u,dudx)where c, f, and s correspond to the , , and terms in Equation 10-2. function [c,f,s] = pdex4pde(x,t,u,DuDx)c = [1; 1];f = [0.024; 0.17] .* DuDx;y = u(1) - u(2);F = exp(5.73*y)-exp(-11.47*y);s = [-F; F];3 Code the initial conditions function. The initial conditions function must be of the formu = icfun(x)The code below represents the initial conditions in the function pdex4ic. function u0 = pdex4ic(x);u0 = [1; 0];4 Code the boundary conditions function. The boundary conditions functions must be of the form[pl,ql,pr,qr] = bcfun(xl,ul,xr,ur,t)The code below evaluates the components p(x,t,u) and q(x,t) (Equation 10-4) of the boundary conditions in the function pdex4bc.function [pl,ql,pr,qr] = pdex4bc(xl,ul,xr,ur,t)pl = [0; ul(2)];ql = [1; 0];pr = [ur(1)-1; 0];qr = [0; 1];5 Select mesh points for the solution. The solution changes rapidly for small t . The program selects the step size in time to resolve this sharp change, but to see this behavior in the plots, output times must be selected accordingly. There are boundary layers in the solution at both ends of [0,1], so mesh points must be placed there to resolve these sharp changes. Often some experimentation is needed to select the mesh that reveals the behavior of the solution.x = [0 0.005 0.01 0.05 0.1 0.2 0.5 0.7 0.9 0.95 0.99 0.995 1];t = [0 0.005 0.01 0.05 0.1 0.5 1 1.5 2];6 Apply the PDE solver. The example calls pdepe with m = 0, the functions pdex4pde, pdex4ic, and pdex4bc, and the mesh defined by x and t at which pdepe is to evaluate the solution. The pdepe function returns the numerical solution in a three-dimensional array sol, wheresol(i,j,k) approximates the kth component of the solution, μk, evaluated at t(i) and x(j).m = 0;sol = pdepe(m,@pdex4pde,@pdex4ic,@pdex4bc,x,t);7 View the results. The surface plots show the behavior of the solution components. u1 = sol(:,:,1); u2 = sol(:,:,2); figure surf(x,t,u1) title('u1(x,t)') xlabel('Distance x') 其输出图形为Distance xu1(x,t)Time tfigure surf(x,t,u2) title('u2(x,t)') xlabel('Distance x') ylabel('Time t')Distance xu2(x,t)Time tAdditional ExamplesThe following additional examples are available. Type edit examplename to view the code and examplename to run the example.。