• MATLAB提供了2种PDE解法，一种是pdepe()函数，他可以求解一般的peds，具有较大的通用性，但只支持命令形式调用；二是pde工具箱，但可以求解特殊PDE问题。
• 非稳态的偏微分方程组是一个比较难解决的问题，也是在热质交换等方面的常常遇到的问题，因此需要一套程序来解决非稳态偏微分方程组的数值
• MATLAB算法-求解微分方程数值解和解析解.ppt
• function [u,x,y] = poisson(f,g,bx0,bxf,by0,byf,D,Mx,My,tol,MaxIter)%solve u_xx + u_yy + g(x,y)u = f(x,y)% over the region D = [x0,xf,y0,yf] = {(x,y) |x0 <= x <= xf, y0 <= y <...

function [u,x,y] = poisson(f,g,bx0,bxf,by0,byf,D,Mx,My,tol,MaxIter)
%solve u_xx + u_yy + g(x,y)u = f(x,y)
% over the region D = [x0,xf,y0,yf] = {(x,y) |x0 <= x <= xf, y0 <= y <= yf}
% with the boundary Conditions:
% u(x0,y) = bx0(y), u(xf,y) = bxf(y)
% u(x,y0) = by0(x), u(x,yf) = byf(x)
% Mx = # of subintervals along x axis
% My = # of subintervals along y axis
% tol : error tolerance
% MaxIter: the maximum # of iterations
x0 = D(1); xf = D(2); y0 = D(3); yf = D(4);
dx = (xf - x0)/Mx; x = x0 + [0:Mx]*dx;
dy = (yf - y0)/My; y = y0 + [0:My]’*dy;
Mx1 = Mx + 1; My1 = My + 1;
%Boundary conditions
for m = 1:My1, u(m,[1 Mx1])=[bx0(y(m)) bxf(y(m))]; end %left/right side
for n = 1:Mx1, u([1 My1],n) = [by0(x(n)); byf(x(n))]; end %bottom/top
%initialize as the average of boundary values
sum_of_bv = sum(sum([u(2:My,[1 Mx1]) u([1 My1],2:Mx)’]));
u(2:My,2:Mx) = sum_of_bv/(2*(Mx + My - 2));
for i = 1:My
for j = 1:Mx
F(i,j) = f(x(j),y(i)); G(i,j) = g(x(j),y(i));
end
end
dx2 = dx*dx; dy2 = dy*dy; dxy2 = 2*(dx2 + dy2);
rx = dx2/dxy2; ry = dy2/dxy2; rxy = rx*dy2;
for itr = 1:MaxIter
for j = 2:Mx
for i = 2:My
u(i,j) = ry*(u(i,j + 1)+u(i,j - 1)) + rx*(u(i + 1,j)+u(i - 1,j))...
+ rxy*(G(i,j)*u(i,j)- F(i,j)); %Eq.(9.1.5a)
end
end
if itr > 1 & max(max(abs(u - u0))) < tol, break; end
u0 = u;
end
以上是possion.m文件，下面给个例子。
%solve_poisson in Example 9.1
f = inline(’0’,’x’,’y’); g = inline(’0’,’x’,’y’);
x0 = 0; xf = 4; Mx = 20; y0 = 0; yf = 4; My = 20;
bx0 = inline(’exp(y) - cos(y)’,’y’); %(E9.1.2a)
bxf = inline(’exp(y)*cos(4) - exp(4)*cos(y)’,’y’); %(E9.1.2b)
by0 = inline(’cos(x) - exp(x)’,’x’); %(E9.1.3a)
byf = inline(’exp(4)*cos(x) - exp(x)*cos(4)’,’x’); %(E9.1.3b)
D = [x0 xf y0 yf]; MaxIter = 500; tol = 1e-4;
[U,x,y] = poisson(f,g,bx0,bxf,by0,byf,D,Mx,My,tol,MaxIter);
clf, mesh(x,y,U), axis([0 4 0 4 -100 100])

展开全文
• 这些由偏微分方程及边界条件、初始条件等组合成的数学模型，只有在十分特殊的条件下才能求得解析解。因此，在很长一段时间内，人们对于这一类问题是无能为力的。随着计算机技术的发展，各种数值方法应运而生，如有限...

【实例简介】
工程中有许多问题可以归结为偏微分方程问题，如弹塑性力学中研究对象(结构、边坡等)内部的应力应变问题、地下水渗流问题等。这些由偏微分方程及边界条件、初始条件等组合成的数学模型，只有在十分特殊的条件下才能求得解析解。因此，在很长一段时间内，人们对于这一类问题是无能为力的。随着计算机技术的发展，各种数值方法应运而生，如有限元法、有限差分法、离散元法、拉格朗日元法等等。利用数值法，可以求得这些问题的数值解。它不是问题的精确解，但可以无限接近精确解。Matlab采用有限元法求解偏微分方程的数值解
偏微吩方程数宽解的 Matlab罢砚
foxdog制作
中输入公式:C1-C2。
3.3定义边界条件
在 Boundary菜单中选择 Boundary mode选项或直接在工具条中单击a按钮,则显示几何模型的外
边界和内边界。如图24所示。
回区
-1;
11l=3,re,命,=|÷正m=:2-1
〔!1〔2〕001〔E
图24定义边界条件
在图中单击边界,则边界线变成黑色,表示它被选中。按住 Shift键,可以连续选择多条边界线段。
选择要定义的边界以后,双击边界或在 Specify Boundary conditions…近项,打开 Boundary Condition
对话框,如图25所示。
吧ndar丁 Condition
瓦[
F IrI dNrs I: ' ll 1(u-I
cue
11
IuI-H
图25 Boundary Condition对话框
对话框中各选项的意义分别为:
Boundary cond. cquation标签显示边界条件方程。
Condition type在下面的两个单选钮中进行选择,确定边界条件类型。
Neumann单选钮选择此项,将边界条件类型确定为 Neuman条仵。
Dirichlet单选钮缺省时选此项。选择此项时,将边界条件类型确定为
Dirichlet边界条件。
● Coelficient value栏在该栏中的对应文丕框中输入边界条件公式卬的系数值。
本例中选择缺省设置。
在 Boundary菜单中近择 Show Edge labels选项和 Show Subdomain labels选项,可以在图中显示
边缘标签和子域标签。
3.4定义PDE类型和PDE系数
偏微分方覆憝宽解的 Matlab寥砚
foxdog制作
在PDE菜单中选择 PdE mode选项,图形将显示为PE模式。选择 Show subdomain labels选项,将
在图中显示子域标签。
在工具条中单击r按钮或在PDE菜单中单击 PDE Specification…选项,可以打开PDE
Specification对话框,如图26所丌。
E日
+改
二1三
图26 PDE Specification对话框
对话框中各选项的意义分别为
Equation标签该标签显示PDE方程表达式。
Type of pD栏在该栏中进行选择,确定PDE问题的类型。
▲ Elliptic单选钮为缺省选项。选择此项,设置为椭圆型PD问题。
Parabolic单选钮选摻此项,设置为抛物线型PDE问题。
▲ Hyperbolic单选钮选择此项,设置为双曲线型PDE问题。
Eigenmodes单选钮选择此项,设置为特征值PE问题。
● Coefficient栏在该栏中设置PDE方程的系数值。
3.5三角形网格剖分
在工具烂中单击△按钮,或在Mesh菜单中选择 Initialize mesh选项,可以进行研究域的三角形
网格初始化。如图27所示。
四17训n小网
上
[吕
-s「∵斗-1.三卩U..甲亠[E了
11
在工具条中单击按钮或在Mesh菜单中选择 Refine mesh选项,可以对初始网格迸行加密。如图
28所示。利用加密后的网格进行计算,可以获得具更高精度的解。
11|」|
偏微吩方程数宽解的 Matlab罢砚
foxdog制作
在Mesh菜单中选择 Jiggle mesh选项,可以对网格进行微调。选择 Display triangle图29微调网
格
Quality,将显示三角形的质量图。如图29所示。图中,将三角形的最好质量定为1,并用红色表
,最差质量为0,用兰色表示,二者之间的三角形质量根据质量得分的大小用红色和兰色的过度色表
小。
在Mesh菜单中选择 Show node labels近项,将在图中显示网格节点的编号,如图30所示。选择
Show Subdomain labels选项,将显示各子域的编号。
其司
C
性
勺
h
12111
图30显示子域编号
在Mesh菜单中选择 Parameters…选项,打开 Mesh parameters刘话框,如图31所示。在该对话框
中进行设置,可以明确与网格剖分有关的参数。
Hesh P
Ah rate
M Jiggle mesh
ggm∈ sh parameters
Jiggle mod
Number of jiggle iterat ons
Refinement method
oK
图31 Mesh Parameters对话框
对话框叩各控件的意义分别为
Maximun edge siκe文本框在该文本框中输入数值,确定最大边缘大小。
Mesh growth rate文本框在该文本框中输入网格生长率。缺刍值为1.3。
Jiggle mesh核选框迒择此项,微调网格。缺省时选择此项。
Jiggle mode下拉式列表栏在其中进行选择,确定微调模式。可选模式有on、 optimize minimum
和 optimize mean等三项。
Number of jiggle iterations文木框在该文本框口输入微调迭代的次数。
Refinement method下拉式列表框在该控件中进行冼择,确定进行网格加密的方法,包括
regular和 longest两个选项。具体内容参见前面函数部分。
偏微吩方程数宽解的 Matlab罢砚
foxdog制作
3.6PDE求解
在工具条中单击=按钮或在 Solve菜单口选择 Solve pDe选项,可以对前面定义的PDE问题进行求
解。本例的解如图32所示。
斗
1·3111
图32问题的缺省图解(色谱图)
在 Solve菜单中单击 Parameters…选项,打开 Solve parameters对话框,如图33所示。在该对话
框中进行设置,可以确定求解方法和参数。
rIlI Ir H■"l
hlurIl'Itri ruri t-=ul lI-'iyle
Ncnlliee l0ler3T
1E4
fuyi'nLrI rur=·urk'l≌.
I'TilinI lill
IIuII_.'Iuu
laurI
FI:h nr r h
J
图33 Solve Parameters对话框
该对话框中各选项的意义分别为
daptive mode核选框选择此项,即选择自适应模式,系统将自适应生成网格并进行求解。
● Maximum number of triangles文本框在其中输入三角形的最大个数。
Maximum number of refinements文本框在其中输入网格加密的最大次数。
Triangle selection method在下面的三个单选钮中进行选择,确定三角形的选择方法。
Worst triangles单选钮为缺省选项。选择此项,根据最坏三角形进行迒择。
▲ Relate tolerance单选钮选此项,根据相对容限进行选择。
▲User- defined function单选钮选择此项,在下面的文本框中输入函数,用该函数选择三角
形
Worst triangle fractions文本框在该文本框中翰入最坏三角形分量。
Refinement method下拉式列表框在其中进行选择,确定进行网格加密的方法,有 regular
和 longest两个选项。
User nonlincar solver核选框选择此项,用非线性求解器进行求解。
偏微吩方程数宽解的 Matlab罢砚
foxdog制作
onl inear tolerance文本框在该文本中输入非线性迭代终止的容限。缺省值为lE-4。
Initial solution文本框在该文本框中输入问题的解的初值。
Jacobian下拉式列表框在该控件中进行选择,确定雅可比矩阵的确定方式。有 fixed、 lumped
和fu11三个选项。
Norm文本框在其口输入范数值,缺省值为Inf(∞)。
3.7解的图形表达
PDE工具项提供了解的多种图形表达方式,上面显示的是彩色图,为缺省显示。单击顽按钮或在
Solve菜单中单击 Parameters…选项,可以打开 Plot selection刈话框,如图34所示。
图34P1 ot selection对话框
F|t/∈
P=:s;
tional
厂 AlmIr,Ad ITE:
leigl H plot I
cor:i I5 T
厂 animato
甲Flls-』 o- SucIm.l
TIn 1.nn
该对话框中各选项的意义分别为:
Plot type控件列该列控件控制图形类型的选择。包括
▲ Color核选框选择此项,生成并显示解的彩色图。缺省时选择此项。
Contour核选框选择此项,生成并显示解的等值线图。如图35所示。
Δ rrows核选框选择此项,生成并显示解的矢量图。如图36所示
De formed mesh核选框选择此项,生成并显示解的变形Ⅸ格图。如图37所示。
▲ Height(3- D plot)核选框选择此项,生成并显示解的三维图。如图38所示。
Animation核选框选择此顶,生成解的系列演示图。
Property空件列该列控件为一组下拉式列表框,定义对解的哪一部分进行图形显小。
User entry控件列为一组文本框,在其中输入用户输入
Plot style控件列该控件列控制前面选择图型的不冋风格。
Plot in x- y grid核选框选择此项,在x-y网格中绘图。
how mesh核选框选择此项,在当前图中显示网格。
Color: u
D
0
E
0
4
0
0.3
"
1-08-0.3-14-0.2C320.40.603112141161
偏微吩方程数宽解的 Matlab罢砚
foxdog制作
图35色谱图叠加等值线图
.2
4
U.2
06
n8-3-14∩.2C12自406811214161.8
图36色谱图叠加矢量图
Y
K
l
4
凵2「2:斗
1叫1
图37变形网格图
E
图38三维高程图
偏微吩方程数宽解的 Matlab罢砚
foxdog制作
Contour plot levels文本框在其中输入等值线的水平数。
Colormap下拉式列表框在该控件中选择绘彩色图的颜色。
Plot solution automatically核选框选择此项,系统自动绘解的图形。
4几种常见的偏微分方程数值求解问题
4.1椭圆型问题
4.1.1单位圆盘的泊松方程
泊松方程是最简单的椭圆型PDE问颎。
该问题的公式为
△U=1
边界上U=0。该问题的精确解为:
U/(x,y)=
y
(1)、用图形用户界面计算
在命令窗囗中输入 pdetool命令,选用 Generic scalar模式
1.单击0 ption菜单,选择 add a grid选项,设置“snap-to-grid”特点。单击+按钮面一个圆,
若该图不是标准的单位圆,双击该圆,打开一对话棰,在其中可以指定圆心的糈确位置和半径的
大小\。
2.通过单击2按钮来设置边界模式。分割的几何边界显示出来,并目外边界指定为缺省设置,即
Dirichlet边界条件,u=0。本例中采月缺省设置,若边界条件不同,可以通过双击边界打开
对话框,在其中入对应的边界条件。
3.单击按丑,定义偏微分方稈,该操作打开一对话框,可以在其中定义PDE系数c,a和f。本
例中,它们均为常数:c-1,f-1,a-0。
4.单击△按扭或选择Mesh菜单中的 Initialize mesh选项,初始化显示三角形Ⅸ格。
5.单击公按扭或在Mesh莫单中选择 Refine mesh选项,改进初始网格并显示新网格。
0.8
0.6
0.4
0.2
0
-0.2
-0.4
-0.6
-0.8
-0.5
0.5
图39初始化网格
偏微吩方程数宽解的 Matlab罢砚
foxdog制作
6.单击=按扭进行求解,Mat1ab可以用图形来表示问题的解。单击按钮,打开 Plot selection
对语框。利月该对话框,可以选择不同类型的解图。
11
[IT
图40问题解的色谱图
7.比较数值解与精确解之间的误差。
为 Plot selecting对话框中的 Color选择 Property弹出式菜单中的 User entry。然后在 User entry
编辑区中输入 Matlab表达式u-(1-x.2-y.2)/4。可以获得解的绝对误差的图形表示。
图41解的绝对误差图形
(二)、使用命令行函数
首先必须创建 MATLAB函数,使二维几何模型参数化。
M文件 circles.m返回单位圆盘边界点的坐标。下面是该文件的内容
nbs=4
if nargin==0
x=nbs;%边界线段个数
retur
eI
0000%参数初值
1111%参数终值
10
【实例截图】
【核心代码】

展开全文
• Matlab偏微分方程工具箱求解方法 这一节我们主要用matlab自带的偏微分方程的工具箱函数求解 一.偏微分方程组的matlab求解语句 ​ 该命令用以求解以下的PDEPDEPDE方程式： c(x,t,u,∂u∂x)∂u∂t=x−m∂(xmf(x,t,...
Matlab的偏微分方程工具箱求解方法

这一节我们主要用matlab自带的偏微分方程的工具箱函数求解

一.偏微分方程组的matlab求解语句
​ 该命令用以求解以下的

P

D

E

PDE

方程式：

c

(

x

,

t

,

u

,

∂

u

∂

x

)

∂

u

∂

t

=

x

−

m

∂

(

x

m

f

(

x

,

t

,

u

,

∂

u

∂

x

)

)

∂

x

+

s

(

x

,

t

,

u

,

∂

u

∂

x

)

c(x,t,u,\frac{\partial u }{\partial x})\frac{\partial u}{\partial t } = x^{-m}\frac{\partial (x^mf(x,t,u,\frac{\partial u}{\partial x}))}{\partial x }+s(x,t,u,\frac{\partial u }{\partial x})

​ 其中：

t

∈

[

t

0

,

t

f

]

,

x

∈

[

a

,

b

]

t \in [t_0,t_f],x \in [a,b]

。偏微分方程的初解：

u

(

x

,

t

0

)

=

v

0

(

x

)

u(x,t_0) = v_0(x)

​ 边界条件为：

p

(

x

,

t

,

u

)

+

q

(

x

,

t

)

f

(

x

,

t

,

u

,

∂

u

∂

x

)

=

0

p(x,t,u) +q(x,t)f(x,t,u,\frac{\partial u}{\partial x}) = 0

​ 下面介绍求解此类方程的函数用法：

s

o

l

=

p

d

e

p

e

(

m

,

p

d

e

p

e

,

i

c

f

u

n

,

b

c

f

u

n

,

x

m

e

s

h

,

t

s

p

a

n

,

o

p

t

i

o

n

s

)

;

sol = pdepe(m,pdepe,icfun,bcfun,xmesh,tspan,options);

​

m

:

m:

对称参数。
​

x

m

e

s

h

:

xmesh:

位置向量,

x

m

e

s

h

=

[

x

0

,

x

1

,

.

.

.

x

N

]

,

x

0

=

a

,

x

N

=

b

xmesh = [x_0,x_1,...x_N],x_0 = a,x_N = b

。
​

t

s

p

a

n

:

tspan:

时间变量

t

t

的向量，

t

s

p

a

n

=

[

t

0

,

t

1

,

.

.

.

t

M

]

,

t

0

=

t

0

,

t

M

=

t

f

tspan = [t_0,t_1,...t_M],t_0 = t_0,t_M = t_f

。
​

p

d

e

f

u

n

:

pdefun:

用户提供的

p

d

e

pde

函数文件。函数格式如下:

[

c

,

f

,

s

]

=

p

d

e

f

u

n

(

x

,

t

,

u

,

d

u

d

x

)

;

[c,f,s] = pdefun(x,t,u,dudx);

​ 也就是说我们要自己设置相应的输出

c

,

f

,

s

c,f,s

，且它们都是行向量。
​

i

c

f

u

n

:

icfun:

求解

u

u

的起始值，格式为

u

=

i

c

f

u

n

(

x

)

u = icfun(x)

。且

u

u

是行向量。
​

b

c

f

u

n

:

bcfun:

提供边界条件函数，格式：

[

p

l

,

q

l

,

p

r

,

q

r

]

=

b

c

f

u

n

(

x

l

,

u

l

,

x

r

,

u

r

,

t

)

;

[pl,ql,pr,qr] = bcfun(xl,ul,xr,ur,t);

​

p

l

,

q

l

:

pl,ql:

左边界

p

p

和

q

q

的行向量。

p

r

,

q

r

:

pr,qr:

右边界

p

p

和

q

q

的行向量。
​

o

p

t

i

o

n

s

:

options:

求解器相关解法参数，见

o

d

e

s

e

t

odeset

。
​

s

o

l

:

sol:

多维向量输出，

s

o

l

(

:

,

:

,

i

)

sol(:,:,i)

为

u

i

u_i

的输出，而

u

i

(

j

,

k

)

=

s

o

l

(

j

,

k

,

i

)

u_i(j,k) = sol(j,k,i)

表示在

t

=

t

s

p

a

n

(

j

)

,

x

=

x

m

e

s

h

(

k

)

t = tspan(j),x = xmesh(k)

时候 的

u

i

u_i

的值。
​ 要获得特定位置和时间的解用以下命令：

[

u

o

u

t

,

d

u

o

u

t

d

x

]

=

p

d

e

v

a

l

(

m

,

x

m

e

s

h

,

u

i

,

x

o

u

t

)

;

[uout,duoutdx] = pdeval(m,xmesh,ui,xout);

​

x

m

e

s

h

:

[

x

0

,

x

1

,

.

.

.

x

N

]

xmesh:[x_0,x_1,...x_N]

​

u

i

:

s

o

l

(

j

,

:

,

i

)

,

ui:sol(j,:,i),

第

i

i

个输出

u

i

u_i

在时间

t

j

t_j

处的解。
​

u

o

u

t

:

uout:

在指定

t

f

t_f

下对应指定位置

x

o

u

t

xout

的值。
​

d

u

o

u

t

d

x

:

duoutdx:

相对应的

d

u

d

x

\frac{du}{dx}

。
二.具体的用法
​ 1.求解以下偏微分方程(解析解为

u

(

x

,

t

)

=

e

−

t

s

i

n

(

π

x

)

u(x,t) = e^{-t}sin(\pi x)

)：

π

2

∂

u

∂

t

=

∂

2

u

∂

x

2

\pi^2\frac{\partial u}{\partial t} = \frac{\partial^2u}{\partial x^2}

​ 其中

x

∈

[

0

,

1

]

x \in[0,1]

,满足以下条件:

u

(

x

,

0

)

=

s

i

n

(

π

x

)

u

(

0

,

t

)

=

0

π

e

−

t

+

∂

u

(

1

,

t

)

∂

x

=

0

u(x,0) = sin(\pi x)\\ u(0,t) = 0 \\ \pi e^{-t} + \frac{\partial u(1,t)}{\partial x} = 0

​

s

o

l

v

e

:

solve:

​ 改写以上偏微分方程到标准形式：

π

2

∂

u

∂

t

=

x

0

∂

∂

x

(

x

0

∂

u

∂

x

)

+

0

\pi^2 \frac{\partial u}{\partial t} = x^0\frac{\partial}{\partial x}(x^0\frac{\partial u}{\partial x}) +0

​ 具体的实现代码如下:
function first
%计算从t:0~3的值
x = linspace(0,1,20);
t = linspace(0,3,60);
subplot(121);
sol = pdepe(0,@firstPdefun,@firstIcfun,@firstBcfun,x,t);
u = surf(x,t,sol(:,:,1));
title('微分方程数值解');
xlabel('x');
ylabel('t');
zlabel('u')
subplot(122);
[X,T] = meshgrid(x,t);
U = exp(-T).*sin(pi*X);
surf(X,T,U);
title('微分方程解析解');
end

%方程段
function [c,f,s] = firstPdefun(x,t,u,dudx)
c = pi^2;
f = dudx;
s = 0;
end
%起始值条件段
function u = firstIcfun(x)
u = sin(pi*x);
end
%边界条件段
function [pl,ql,pr,qr] = firstBcfun(xl,ul,xr,ur,t);
pl = ul;
ql = 0;
pr = pi*exp(-t);
qr = 1;
end

对比一下发现几乎和解析解一摸一样！
2.求解以下偏微分方程的数值解：

∂

u

1

∂

t

=

0.024

∂

2

u

1

∂

x

2

−

F

(

u

1

−

u

2

)

∂

u

2

∂

t

=

0.170

∂

2

u

2

∂

x

2

+

F

(

u

1

−

u

2

)

F

(

u

1

−

u

2

)

=

e

5.73

(

u

1

−

u

2

)

−

e

−

11.46

(

u

1

−

u

2

)

\frac{\partial u_1}{\partial t} = 0.024\frac{\partial^2u_1}{\partial x^2} - F(u_1-u_2)\\ \frac{\partial u_2}{\partial t} = 0.170\frac{\partial^2u_2}{\partial x^2} + F(u_1-u_2)\\ F(u_1-u_2) = e^{5.73(u_1-u_2)} - e^{-11.46(u_1-u_2)}

初值条件:

u

1

(

x

,

0

)

=

1

u

2

(

x

,

0

)

=

0

u_1(x,0) = 1\\ u_2(x,0) = 0

边值条件:

∂

u

1

(

0

,

t

)

∂

x

=

0

u

2

(

0

,

t

)

=

0

u

1

(

1

,

t

)

=

1

∂

u

2

(

1

,

t

)

∂

x

=

0

\frac{\partial u_1(0,t)}{\partial x} = 0\\ u_2(0,t) = 0\\ u_1(1,t) = 1\\ \frac{\partial u_2(1,t)}{\partial x} = 0\\

S

o

l

v

e

:

Solve:

​ 化简上面的偏微分方程为标准形式：

(

1

1

)

.

∗

∂

∂

t

(

u

1

u

2

)

=

∂

∂

x

(

0.024

∂

u

1

∂

x

0.170

∂

u

2

∂

x

)

+

(

−

F

(

u

1

−

u

2

)

F

(

u

1

−

u

2

)

)

\begin{pmatrix} 1\\1 \end{pmatrix}_.*\frac{\partial}{\partial t}\begin{pmatrix} u_1\\u_2 \end{pmatrix} = \frac{\partial}{\partial x}\begin{pmatrix} 0.024\frac{\partial u_1}{\partial x}\\0.170\frac{\partial u_2}{\partial x} \end{pmatrix} +\begin{pmatrix} -F(u_1-u_2)\\F(u_1-u_2) \end{pmatrix}

​ 化简左边界条件也有：

(

0

u

2

)

+

(

1

0

)

.

∗

(

0.024

∂

u

1

∂

x

0.170

∂

u

2

∂

x

)

=

(

0

0

)

\begin{pmatrix} 0\\u_2 \end{pmatrix}+\begin{pmatrix} 1\\0 \end{pmatrix}_.* \begin{pmatrix} 0.024\frac{\partial u_1}{\partial x}\\0.170\frac{\partial u_2}{\partial x} \end{pmatrix}= \begin{pmatrix} 0\\0 \end{pmatrix}

​ 化简右边界条件有：

(

u

1

−

1

0

)

+

(

0

1

)

.

∗

(

0.024

∂

u

1

∂

x

0.170

∂

u

2

∂

x

)

=

(

0

0

)

\begin{pmatrix} u_1-1\\0 \end{pmatrix}+\begin{pmatrix} 0\\1 \end{pmatrix}_.* \begin{pmatrix} 0.024\frac{\partial u_1}{\partial x}\\0.170\frac{\partial u_2}{\partial x} \end{pmatrix}= \begin{pmatrix} 0\\0 \end{pmatrix}

​ 最后附上整个代码：
function second
xmesh  = linspace(0,1,20);
tspan = linspace(0,3,60);
sol = pdepe(0,@secondPdefun,@secondPdein,@secondPdebc,xmesh,tspan);
subplot(121);
surf(xmesh,tspan,sol(:,:,1));
xlabel('x');
ylabel('t');
zlabel('u(1)');
title('u(1)-x-t图像');
subplot(122);
surf(xmesh,tspan,sol(:,:,2));
xlabel('x');
ylabel('t');
zlabel('u(2)');
title('u(2)-x-t图像');
end

function [c,f,s] = secondPdefun(x,t,u,dudx)
c = [1 1]';
f = [0.024*dudx(1) 0.170*dudx(2)]';
y = u(1) - u(2);
stemp = exp(5.73*y) - exp(-11.46*y);
s = [-stemp stemp]';
end

function up = secondPdein(x,t,u,dudx)
up = [1 0]';
end

function [pl ql pr qr] = secondPdebc(xl,ul,xr,ur,t)
pl = [0 ul(2)]';
ql = [1 0]';
pr = [ur(1)-1 0]';
qr = [0  1]';
end


最后的结果是：

展开全文
• 文章目录微分方程概述及matlab代码（求解析解微分方程概述引例：导弹射击问题微分方程基本概念Matlab微分方程解析解关于导弹射击问题的matlab求解关于Matlab解析解的几点注意 微分方程概述 当我们描述实际...
微分方程概述及matlab代码（求解析解）
【清风数学建模教程笔记】

文章目录
微分方程概述及matlab代码（求解析解）微分方程概述引例：导弹射击问题微分方程基本概念Matlab求微分方程的解析解关于导弹射击问题的matlab求解关于Matlab求解析解的几点注意

微分方程概述
当我们描述实际对象的某些特性随时间（空间）而演变的过程，分析它的变化规律，预测它的未来形态、研究它的控制手段时。通常要建立对象的动态模型。
在许多实际问题中，当直接导出变量之间的函数关系较为困难，但导出包含未知函数的导数或微分的关系式比较容易时，可以用建立微分方程模型的方法来研究该问题。
分类：常微分方程（未知函数是一个一元函数）和偏微分方程（未知函数是一个多元函数）国赛趋势：越来越重视微分方程、物理领域等解析解和数值解
解析解(analytical solution)是严格按照公式逻辑推导得到的，具有基本的函数形式。给出任意的自变量就可以求出其因变量，也就是问题的解，他人可以利用这些公式计算各自的问题，具有广泛适用性；数值解(numerical solution)是采用某种计算方法，在特定的条件下得到的一个近似数值结果，如有限元法，数值逼近法，插值法等等得到的解。别人只能利用数值计算的结果，而不能随意给出自变量并求出计算值。 微分方程模型的应用十分广泛：例如人口预测模型、捕食者猎物模型、种群相互竞争\依存模型、传染病模型等
关于传染病模型的学习资源
上海财经大学传染病模型讲解上海市现代应用数学重点实验室对新冠肺炎的一些研究报告(emo 看不懂)
引例：导弹射击问题

微分方程基本概念
微分方程：含导数或微分的方程称为微分方程，其一般形式为

f

(

x

,

y

1

,

y

2

,

.

.

.

,

y

n

)

=

0

f(x,y^1,y^2,...,y^{n})=0

微分方程的阶数：微分方程中所含导数或微分的最高阶数称为微分方程的阶数，例如下式是三阶微分方程。

y

′

′

′

+

2

y

′

′

−

2

x

=

0

y'''+2y''-2x = 0

微分方程的解：使得微分方程成立的函数称为微分方程的解。  微分方程的通解和特解：不含任意常数的解称为微分方程的特解；若解中所含的相互独立的任意常数的个数与微分方程的阶数相等，则称此解为微分方程的通解。例如：

y

′

−

e

x

=

0

y'-e^x=0

特解为

e

x

e^x

通解为

e

x

+

C

e^x+C

初值条件：能够确定通解中任意常数的条件称为微分方程的初值条件。  如何建立微分方程？
按照专业知识、规律直接列方程。套用现有的微分方程模型。（根据具体问题改进）
Matlab求微分方程的解析解
解析解：即给出解的具体表达式，若给定初始条件，则能求出含具体常数的解。
以下给出几个例子。
例1：

y

−

d

y

d

x

=

2

x

y - \frac{dy}{dx}=2x

clear;clc
dsolve('y-Dy=2*x','x')
%此处要指定自变量为x
% 输出 2*x + C1*exp(x) + 2  (这里的C1表示任意常数，有时候也会出现C2 C3等)
% 如果不指定自变量的话，会默认自变量为t，x会看成一个常数

% 下面这种写法是新版的matlab推荐的方式,注意：最新版本的matlab会逐渐淘汰上面那种写法
syms y(x)
eqn = (y - diff(y,x) == 2*x);    % 注意原来方程中的“=”一定要改成“==”
dsolve(eqn)

当微分方程中还有其他的未知参数时

y

−

d

y

d

x

=

a

x

y - \frac{dy}{dx}=ax

% 方法1
dsolve('y-Dy=a*x','x')  % a是一个未知的参数
% 方法2
syms y(x) a
eqn = (y - diff(y,x) == a*x);
dsolve(eqn)

例2：

y

−

y

′

=

2

x

且

y

(

0

)

=

3

y-y' =2x且y(0)=3

% 方法1
dsolve('y-Dy=2*x','y(0)=3','x')
% 输出：2*x + exp(x) + 2
% 方法2
syms y(x)
eqn = (y - diff(y,x) == 2*x);
%初值条件
cond = (y(0) == 3);
dsolve(eqn,cond)
% 2*x + exp(x) + 2

例3：

y

′

′

+

4

y

′

+

29

y

=

0

且

y

(

0

)

=

0

,

y

′

(

0

)

=

15

y'' +4y'+29y=0且y(0)=0,y'(0)=15

% 方法1
dsolve('D2y+4*Dy+29*y=0','y(0)=0,Dy(0)=15','x')
% 3*sin(5*x)*exp(-2*x)
% 方法2
syms y(x)
eqn = (diff(y,x,2) + 4 *diff(y,x) + 29*y  == 0);
Dy = diff(y,x); % 定义变量Dy为y的一阶导数
cond = [(y(0) == 0) ,(Dy(0) ==15)] ;
% 有两个条件，可以写到一个向量中保存
dsolve(eqn,cond)
% 3*sin(5*x)*exp(-2*x)

例4：
微分方程组求解，具体见以下代码。
% 方法1
[x,y,z] = dsolve('Dx=2*x-3*y+3*z+t','Dy=4*x-5*y+3*z+t','Dz=4*x-4*y+2*z+t','t')

% 方法2
syms x(t) y(t) z(t)
eqn1 = (diff(x,t)  == 2*x-3*y+3*z+t);
eqn2 = (diff(y,t)  == 4*x-5*y+3*z+t);
eqn3 = (diff(z,t)  == 4*x-4*y+2*z+t);
eqns = [eqn1 eqn2 eqn3];
[x,y,z] = dsolve(eqns)
% x = exp(2*t)*(C2- (exp(-2*t)*(2*t + 1))/4) + C3*exp(-t)
% y = exp(2*t)*(C2 - (exp(-2*t)*(2*t + 1))/4) + C3*exp(-t) + C4*exp(-2*t)
% z = exp(2*t)*(C2 - (exp(-2*t)*(2*t + 1))/4) + C4*exp(-2*t)

关于mupad工具包以及不手打公式的工具等
mupad  % 最新版本matlab可能会报错，将计算结果复制到里面，使结果可读。
% 如果新版matlab用不了mupad的话，可以使用实时脚本（【新建】-【实时脚本】）
simplify(y)  % simplify函数可以简化表达式
latex(y)
% 转换成latex代码，复制到Axmath或者word自带的公式编辑器
% 如果太过于复杂的话可能会报错，可以自己测试

关于导弹射击问题的matlab求解
% 假设 v=100
[x,y] = dsolve('Dx = 3*100*(20+sqrt(2)/2*100*t-x)/sqrt((20+sqrt(2)/2*100*t-x)^2+(sqrt(2)/2*100*t-y)^2)','Dy = 3*100*(sqrt(2)/2*100*t-y)/sqrt((20+sqrt(2)/2*100*t-x)^2+(sqrt(2)/2*100*t-y)^2)','x(0)=0,y(0)=0','t')
% 警告: Explicit solution could not be found. 不能找到显示解

此题无法使用求解析解的方法求解。
关于Matlab求解析解的几点注意
自变量需要指定，当不指定时默认为t初始条件不给时，算出来的是通解可以用向量的方式保存微分方程组并求解如果微分方程形式较为复杂，那么很大可能是得不到解析解的，则需要求数值解
展开全文
• 这部分主要讨论如何用MATLAB实现对偏微分方程的数值仿真求解．MATLAB偏微分方程工具箱(PDE Toolbox)的出现，为偏微分方程的求解以及定性研究提供了捷径．主要步骤为： 2.1 用偏微分方程工具箱求解微分方程 直接...
• 偏微分方程的数值(二): 一维状态空间的偏微分方程MATLAB 解法 偏微分方程的数值(三): 化工应用实例 ----------触煤反应装置内温度及转换率的分布 偏微分方程的数值(四): 化工应用————扩散系统之浓度...
• 二阶椭圆偏微分方程实例求解(附matlab代码).docx 《微分方程数值解法》期中作业实验报告二阶椭圆偏微分方程第一边值问题姓名：学号：班级：2013年11月19日1二阶椭圆偏微分方程第一边值问题摘要对于二阶椭圆偏微分...
• 零基础使用 MATLAB 求解偏微分方程（建议收藏） ...在所有的偏微分方程中，百分之九十九都是没有解析解的。没有解析解怎么办，我们只能通过有限元或者有限差分等方法，求解偏微分方程数值解。如果您有一些代码基础，
• 偏微分方程的数值(二): 一维状态空间的偏微分方程MATLAB 解法 偏微分方程的数值(三): 化工应用实例 ----------触煤反应装置内温度及转换率的分布 偏微分方程的数值(四): 化工应用————扩散系统之浓度...
• 之一：【MATLAB】逐步搜索法、二分法、比例求根法、牛顿法、弦截法求方程的根之二：【MATLAB】...Nicholson隐式法(抛物型偏微分方程) 上图疑似有误，应为 Crank-Nicholson 隐式法边值为u显式法 % 显式法(抛物型偏微...
• matlab微分方程实验报告4《matlab与数学实验》实验报告实验序号： 实验四 日期： 2015年 5 月 25 日班级 132132002姓名 彭婉婷学号 1321320057实验名称 求微分方程 问题背景描述实际应用问题通过数学建模...
• 第 一 卷 第 期 河 北理 工大学学报 自然科学 版年 月 ‘ 【 立 ·‘爹编号 扳 抖 刁 一 阶偏微分方程 的 和 解法 比较 谷建涛 , 终 玉 霞 , 付景红 河北理工大学 理学院 河北 脚 肠 沟 关健询 一阶偏微分 方程 摘 ...
• 针对一阶线性双曲型偏微分方程,要求其Cauchy问题的解析解,提出特征线方法 .特征线方法的基本思想是将偏微分方程的Cauchy问题转化为常微分方程的相应问题,通过解常微分方程进而得到原来偏微分方程问题的解.通过对特征...
• 微分方程Matlab绘图微分方程Matlab绘图2009.07Edit by niubenEdit by niubenContentsContents1 Matlab绘图问题1 Matlab绘图问题2 常微分方程求解问题2 常微分方程求解问题3 微分方程 matlab 求解3 微分...
• ## matlab求解微分方程

千次阅读 2020-03-30 00:51:14
matlab求解微分方程 1.解析解 dsolve 函数 dsolve函数用于求常微分方程组的精确解，也称为常微分方程的符号解。如果没有初始条件或边界条件，则求出通解；如果有，则求出特解。 1)函数格式 Y = dsolve(‘eq1,eq2,…...
• (1)形如 du/dt + a*du/dx = v*du2/dx2(二次偏导数) 的一阶偏微分方程可用pde函数求解如下例：0编程如下：function pdexvvm = 0;x = linspace(0,1,101);t = linspace(0,0.00006,501);sol = pdepe(m,@pdexvvpde,@...
• 《基于MATLAB偏微分方程差分解法》由会员分享，可在线阅读，更多相关《基于MATLAB偏微分方程差分解法(12页珍藏版)》请在人人文库网上搜索。1、基于MATLAB偏微分方程差分解法学院：核工程与地球物理学院专业：...
• 本篇将介绍用matlab求解常微分方程的数值解和解析解，并非是一种完整的模型，仅仅是一些算法。由于数学原理过于复杂，故不探究背后的数学原理，仅将matlab求解的相关函数加以记录。所有代码均可跑通。 1.Matlab求常...
• 1.求解拉普拉斯方程的狄利克雷法 求解在区域R = {(x,y): 0≤x≤a,0≤y≤b}内的 uxx(x,y) + uyy(x,y) = 0 的近似，而且满足条件 u(x,0) = f1(x),u(x,b) = f2(x), 其中0≤x≤a 且u(0,y) = f3(y),u(a,y) = f4(y),...
• 是保守双曲偏微分方程 (PDE) 的非线性系统。 该系统存在黎曼问题的精确，包括冲击波和中心稀疏波的不同组合。 所考虑的经典问题是所谓的溃坝问题，其中水流最初处于静止状态（零速度），水深 h 具有阶梯不连续性，...
• matlab求微分方程的.doc matlab求微分方程的解一、问题...另一方面，能够求解的微分方程也是十分有限的，特别是高阶方程和偏微分方程(组)．这就要求我们必须研究微分方程(组)的解法，既要研究微分方程(组)的解析...
• 偏微分方程，再加上边界条件、初始条件构成的数学模型，只有在很特殊情况下才可求得解析解。随着计算机技术的发展，采用数值计算方法，可以得到其数值解。下面的几个简单例子，将为大家介绍如何利用Matlab中的PDE...
• 偏微分的第4次实验，主要内容为二维抛物型方程的ADI格式求解。不足之处望读者多多指正。 实验内容 使用ADI格式求解以下问题： ∂u∂t−(∂2u∂x2+∂2u∂y2)=−32e12(x+y)−t,0,y,0⩽1\frac{\partial u}{\partial t}-...
• Matlab学习——求解微分方程(组)发布时间：2018-05-29 17:18,浏览次数：738, 标签：Matlab介绍：1.在 Matlab 中，用大写字母 D 表示导数，Dy 表示 y 关于自变量的一阶导数，D2y 表示 y 关于自变量的二阶导数，...
• 1、偏微分方程的,一、离散化的概念,油藏是非均质的，岩石和流体性质伴随时间常常是发生变化的，建立的偏微分方程一般是非线性的，求解偏微分方程解析解比较困难，常用数值求解。 目前工程上应用的离散化方法有：...
• ## matlab解常微分方程

千次阅读 2021-03-09 20:19:52
微分方程ordinary differential equation的缩写，此种表述方式常见于编程，如MATLAB中Simulink求解器solver已能提供了7种微分方程求解方法：ode45(Dormand-Prince)，ode23(Bogacki-Shampine)，ode113(Adams)，ode...

...

matlab 订阅