PDE 有限差分法

u(2:end-1,3-i)=u(2:end1,i)+a*dt./(r.^2).*( r.^2.*(u(3:end,i)+u(1:end-2,i)-2*u(2:end1,i))/dr^2 + r.*(u(3:end,i) - u(1:end-2,i))./(2*dr) - u(2:end1,i) ) + b./r*dt;
i=3-i; t=t+dt; if (mod(round(t/dt),1000)==0)
plot(r0,u(:,i)); title(['time=',num2str(t),' s']); pause(0.1); end end
参考书：顾昌鑫，计算物理学.上海：复旦大学出版社.2010
+
������������−1,������+1
−
2������������−1,������ ℎ2
+
������������−1,������−1]
+
������������,������
这种差分方法对任意离散方法微扰收敛。
思考：为什么 Von Neumann 方法需要求解线性方程组？你能写一个 Matlab 程序实现它吗？
从而我们需要采用新的差分方法—三点式差分：
������������ ������������
=
������������+1,������
−
(������������,������+1 ������
+
������������,������−1)/2
思考：构造相同的微绕，请自己证明收敛条件为|������������| < ������，从而在恰当离散下可以收敛
������������ ������������ |������=0
=
������(������)
{������|������=0 = ������(������) ������|������=������ = ������(������)
(1)三层显示差分：
������������+1,������
+
������
������������,������+1 − ������������,������−1 2ℎ
=
0
在初值条件中引入微扰项ϵk(������) 从而在������|������=0 = ������(������) + ϵk(������)下求解原方程，假设ϵk(������) = exp(������������������������)，在加入微扰后的解为 ���̃���(������, ������)，记v(x, t) = ���̃���(������, ������) − ������(������, ������),则v(x, t)满足的差分方程为：
������������
分都将会导致发散（计算误差的迭代）。
（发散的证明：如果������������采用中心差分，������������采用向前差分，则差分方程为：
������������
������������
������������+1,������ − ������
������������,������
=
������2������2
[������������,������+1
−
2������������,������ ℎ2
+
������������,������−1
+
������������,������ ]
+
2������������,������
−
������������−1,������
������������,0 = ������������
������������ ������������
−
������)
+
������ ������
������|������=0 = 0
{
������|������=0 = 0 ������|������=������ = 0
%% parameter
rmin=0;
rmax=0.01; dr=rmax/100; dt=0.0001; tmax=1000; u0=0; ub=0;
a=1e-6; b=1e-9;
%a=mu/rho %b=Delta/rho
%% init r=(rmin:dr:rmax)'; r0=r; Nr=length(r); u=zeros(Nr,2); u(:,1)=u0; u(1,:)=ub; u(end,:)=ub;
%% calc r=r(2:end-1); t=0; i=1; while (t<=tmax)
������������,1
− ������
������������,0
Leabharlann Baidu
=
������������
{
������0,������ = ������������, ������������,������ = ������������
微扰收敛条件|aτ| < 1
h
(2)Von Neumann 差分：
取 a=1，
(4 一维二阶 PDE 热传导 显示差分、隐式差分、六点隐式差分
5 Poission 方程求解 …)
6 例子：CUPT2015 4 Liquid Film Motor
Matlab 实现：
∂u ∂t
=
������ ������2
(������2
������2������ ������������2
+
������
������������+1,������ −
{
������
������������,������ + ������ ������|������=0 =
������������,������+1 − ������������,������−1
2ℎ exp(������������������������)
������|������=0 = ������1(������) ������|������=������ = ������2(������)
求解步骤：
(1)离散化（恰当选取 h,τ,保证求解结果不发散）
(2)写出差分
(3)写出差分方程
(4)计算（直接由前一层计算后一层）
几种差分方法：向前差分、向后差分、中心差分
=
0
������|������=0 = 0 ������|������=������ = 0
求解得到： 其中
vk,j = ������������exp(������������������������)
aτ λ = 1 − i h sin(������ℎ)
从而λ模长大于 1，从而发散。） 思考：如果两个都采用向前差分如何构造微扰证明发散？
������������采用中心差分
������������
������������ (������������)������,������
=
������������,������+1 − ������������,������−1 2ℎ
������������采用这三种都不行，因为：中心差分将依赖前两层，而初值仅仅给出一层；向前或向后差
������������+1,������ − 2������������,������ + ������������−1,������ ������2
=
1 2
[������������+1,������+1
−
2������������+1,������ ℎ2
+
������������+1,������−1
PDE 有限差分法讲义
2015/11/28 刘语涵
1 偏微分方程简介（方程+初值条件+边界条件） 2 一维一阶 PDE 例子：
������������ ������������������
{
������������ + ������������ = 0 ������|������=0 = ������(������)
������
思考：如何解第二类边界条件？
3 一维二阶 PDE 弦振动
������2������ ������2������2������ ������������2 + ������������2 = ������(������, ������)
������|������=0
=
������(������)