解偏微分方程(研究生课程)

在右侧热量以恒定速率降低到 周围空气中，其它边界独立。 即： ⑴ u=100 左侧（dirichilet条件） ⑵ u’=-10 右侧（Neumann条件） ⑶ u’=0 其它边界（Neumann条件）
Rcdao.m
wk.baidu.com
差分法解热传导方程
一维热传导方程是： ut a2uxx
由导数差分公式有： ut
3.选定方程的类型 Elliptic 椭圆型 f=0 Parabolic 抛物线 Hyperbolic 双曲线 Eigenmodes 特征值 4.将网格初始化和精细化 5.求解方程 6.将结果画图 7.放大视图
（抛物线型）受热金属块的热传导
一块受热的有矩形裂纹的金属块， 金属块的左侧被加热到100摄氏度，
波动方程可以写成差分 形式：
u( x, t

t )
2u( x, t) （t）2

u( x, t

t )

a2
u( x

x, t)

2u( x, t) (x)2

u( x

x,
t)
令: x

ix, t

jt, c

(t )2 (x)2
a2,
i, j 0,1,2,...n 1
u( x, t

t )

u( x,
t)

t (x)2
a2[u( x

x,
t)

2u( x, t)

u( x

x, t)]
令: x

ix, t

jt, r

t (x)2
a2,
i, j 0,1,2,...n 1
ui,j+1
写成足标形式
ui, ji ui, j rui1, j 2rui, j rui1, j (1 2r )ui, j r(ui1, j ui1, j )
% -div(grad(u))=1
U(x, y) 1 x2 y2 4
% on the unit disk with u=0 on the boundary. % Compare with exact solution. pause % Strike any key to continue.
%差分法解热传导方程 cfrcd.m
x=0:20;a2=10;r=a2*0.01;
u=zeros(21,25); %预设矩阵以存放求得的解
u(10:11,1)=1; %初始条件
for j=1:25
%求解及作图
u(2:20,j+1)=(1-2*r)*u(2:20,j)+r*(u(1:19,j)+u(3:21,j));
rui-1,j (1-2r)ui,j
rui+1,j
差分法解热传导方程
有限长细杆的热传导问 题
ut u(0,
a 2uxx t) 0,
u(l
,
t
)

0
u( x, t 0) ( x)
参数：l 20, t 25, a2 10,且

(
x)

1
0
(10 x 11) ( x 10, x 11)
h=plot(x,u(:,1),'linewidth',3);
axis([0,1,-0.05,0.05]);
set(h,'EraseMode','xor','MarkerSize',18)
for k=2:N
set(h,'XData',x,'YData',u(:,2));
drawnow;
u(2:419,3)=2*u(2:419,2)-u(2:419,1)+c*(u(3:420,2)-2*u(2:419,2)+u(1:418,2));
dt=0.0005;c=dt*dt/dx/dx;
x=linspace(0,1,420)';
u(1:420,1)=0;
u(181:240,1)=0.05*sin(pi*x(181:240)*7);
u(2:419,2)=u(2:419,1)+c/2*(u(3:420,1)-2*u(2:419,1)+u(1:418,1));
引言：偏微分方程的主要类型
椭圆型

(
2 x 2

2 y 2

2 z 2
)( x,
y, z)

S( x,
y, z)
抛物型
[ t

a
2
(
2 x 2

2 y 2

2 z 2
)] (
x,
y,
z,
t)

S(
x,
y,
z,
t)
双曲型 2
[t 2

a
2
(
2 x 2

2 y 2
• MATLAB采用有限元法求解偏微分方程的数值解。
偏微分方程工具箱
《数学物理方程的 MATLAB解法与可视化》 彭芳麟 清华大学出版社
详细地介绍了MATLAB 的偏微分方程工具箱与解 偏微分方程指令，还介绍 了差分方法和有限元方法。 对学习数值计算或计算物 理课程而言，这也是很实 用的参考教材。
u(2:419,1)=u(2:419,2);
u(2:419,2)=u(2:419,3);
end
（抛物线型）矩形区域的有源的热传导
在矩形的四边,温度为零,在区域的中央,半径为 0.4的圆内,有个恒定的热源,其值为1.
（双曲线型）方形薄膜的横向振动的波动方程
方形薄膜左侧和右侧固定 （u=0）,前后两端自由（u’=0） 初始条件：t=0
u(0)=atan(cos(pi/2*x)) dudt(0)=3*sin(pi*x).*exp(sin(pi/2*y))

2 z 2
)]( x,
y,
z, t )

S( x,
y, z,t)
拉普拉斯方程 热传导方程、扩散方程 波动方程
偏微分方程工具箱（PDETOOL）
（椭圆型）单位圆盘的泊松方程
在单位圆内求解泊松方程
-ΔU=1
在单位圆的边界上U=0。 % Solve Poisson's equation
该问题的精确解为
偏微分方程的主要类型222椭圆型???????xyz?sxyz222?x?y?z222?2????a???xyzt?sxyzt抛物型222?t?x?y?z2222?2????a???xyzt?sxyzt2222双曲型?t?x?y?z拉普拉斯方程热传导方程扩散方程波动方程偏微分方程工具箱pdetool椭圆型单位圆盘的泊松方程在单位圆内求解泊松方程u1在单位圆的边界上u0
（椭圆型）矩形区域的拉普拉斯方程
任意选取定解问题中参数的值,例 如取μ=1,a=1,b=1 这个问题的解析解如下
（椭圆型）矩形区域的拉普拉斯方程
解方程的操作步骤：
1.画方程求解的区域 2.设定边界条件dirichlet 边界1：sin(3*pi.*x).*cos(pi.*x) 边界2：sin(3*pi.*y)
写成足标形式
ui, ji ui, j rui1, j 2rui, j rui1, j c(ui1, j ui1, j ) 2(1 c)ui, j ui, j1
2-2c
差分法解弦振动方程
%两端固定的弦振动 clear
两端固定的均匀弦的振动
N=4010;dx=0.0024;
数学物理方程的Matlab解法
数学物理方程的Matlab解法
• 工程中许多问题可以归结为偏微分方程问题，这 些由偏微分方程及边界条件、初始条件等构成的 数学模型，只有在十分特殊的条件下才能求得解 析解。
• 随着计算机技术的发展，各种数值方法应运而生， 如有限元法，有限差分法、拉格朗日元法等。利 用数值法，可以求得这些问题的数值解。它不是 问题的精确解，但可以无限接近精确解。
plot(x,u(:,j));axis([0 21 0 1]);pause(0.1)
end
meshz(u)
差分法解弦振动方程
一维弦振动方程： utt a2uxx
初始条件：u(
x,0)

x 1.5

(0 1.5 x
x
3 / 5) (3/ 5
x

1)
ut ( x,0) 0

u( x, t

t) t
u( x, t)
uxx

u( x

x, t)

2u( x, t) (x)2
u( x

x, t)
热传导方程可以写成差 分形式：
u(
x, t

t ) t

u( x,
t)

a2
u( x

x, t)

2u( x, t) (x)2

u( x

x, t)
即：