搜索
PDE有限差分法讲义2015/11/28刘语涵1偏微分方程简介（方程+初值条件+边界条件）2一维一阶PDE例子：{𝜕𝑢𝜕𝑡+𝑎𝜕𝑢𝜕𝑥=0𝑢|𝑡=0=𝜙(𝑥)𝑢|𝑥=0=𝑓1(𝑡)𝑢|𝑥=𝑅=𝑓2(𝑡)求解步骤：(1)离散化（恰当选取h,τ,保证求解结果不发散）(2)写出差分(3)写出差分方程(4)计算（直接由前一层计算后一层）几种差分方法：向前差分、向后差分、中心差分𝜕𝑢𝜕𝑥采用中心差分(𝜕𝑢𝜕𝑥)𝑘,𝑗=𝑢𝑘,𝑗+1−𝑢𝑘,𝑗−12ℎ𝜕𝑢𝜕𝑡采用这三种都不行，因为：中心差分将依赖前两层，而初值仅仅给出一层；向前或向后差分都将会导致发散（计算误差的迭代）。（发散的证明：如果𝜕𝑢𝜕𝑥采用中心差分，𝜕𝑢𝜕𝑡采用向前差分，则差分方程为：𝑢𝑘+1,𝑗−𝑢𝑘,𝑗𝜏+𝑎𝑢𝑘,𝑗+1−𝑢𝑘,𝑗−12ℎ=0在初值条件中引入微扰项ϵk(𝑥)从而在𝑢|𝑡=0=𝜙(𝑥)+ϵk(𝑥)下求解原方程，假设ϵk(𝑥)=exp(𝑖𝑎𝑘𝑗)，在加入微扰后的解为𝑢̃(𝑥,𝑡)，记v(x,t)=𝑢̃(𝑥,𝑡)−𝑢(𝑥,𝑡),则v(x,t)满足的差分方程为：{𝑣𝑘+1,𝑗−𝑣𝑘,𝑗𝜏+𝑎𝑣𝑘,𝑗+1−𝑣𝑘,𝑗−12ℎ=0𝑣|𝑡=0=exp(𝑖𝑎𝑘𝑗)𝑣|𝑥=0=0𝑢|𝑥=𝑅=0求解得到：vk,j=𝜆𝑗exp(𝑖𝑎𝑘𝑗)其中λ=1−iaτhsin(𝑎ℎ)从而λ模长大于1，从而发散。）思考：如果两个都采用向前差分如何构造微扰证明发散？从而我们需要采用新的差分方法—三点式差分：𝜕𝑢𝜕𝑡=𝑢𝑘+1,𝑗−(𝑢𝑘,𝑗+1+𝑢𝑘,𝑗−1)/2𝜏思考：构造相同的微绕，请自己证明收敛条件为|𝐚𝛕𝐡|𝟏，从而在恰当离散下可以收敛思考：如何解第二类边界条件？3一维二阶PDE弦振动{𝜕2𝑢𝜕𝑡2+𝑎2𝜕2𝑢𝜕𝑥2=𝑓(𝑥,𝑡)𝑢|𝑡=0=𝜙(𝑥)𝜕𝑢𝜕𝑡|𝑡=0=𝜓(𝑥)𝑢|𝑥=0=𝑔(𝑡)𝑢|𝑥=𝑅=𝑚(𝑡)(1)三层显示差分：{𝑢𝑘+1,𝑗=𝑎2𝜏2[𝑢𝑘,𝑗+1−2𝑢𝑘,𝑗+𝑢𝑘,𝑗−1ℎ2+𝑓𝑘,𝑗]+2𝑢𝑘,𝑗−𝑢𝑘−1,𝑗𝑢𝑗,0=𝜙𝑗𝑢𝑗,1−𝑢𝑗,0𝜏=𝜓𝑗𝑢0,𝑘=𝑔𝑘,𝑢𝑁,𝑘=𝑚𝑘微扰收敛条件|aτh|1(2)VonNeumann差分：取a=1，𝑢𝑘+1,𝑗−2𝑢𝑘,𝑗+𝑢𝑘−1,𝑗𝜏2=12[𝑢𝑘+1,𝑗+1−2𝑢𝑘+1,𝑗+𝑢𝑘+1,𝑗−1ℎ2+𝑢𝑘−1,𝑗+1−2𝑢𝑘−1,𝑗+𝑢𝑘−1,𝑗−1ℎ2]+𝑓𝑘,𝑗这种差分方法对任意离散方法微扰收敛。思考：为什么VonNeumann方法需要求解线性方程组？你能写一个Matlab程序实现它吗？(4一维二阶PDE热传导显示差分、隐式差分、六点隐式差分5Poission方程求解…)6例子：CUPT20154LiquidFilmMotor{∂u∂t=𝑎𝑟2(𝑟2𝜕2𝑢𝜕𝑟2+𝑟𝜕𝑢𝜕𝑟−𝑢)+𝑏𝑟𝑢|𝑡=0=0𝑢|𝑥=0=0𝑢|𝑥=𝑅=0Matlab实现：%%parameterrmin=0;rmax=0.01;dr=rmax/100;dt=0.0001;tmax=1000;u0=0;ub=0;a=1e-6;%a=mu/rhob=1e-9;%b=Delta/rho%%initr=(rmin:dr:rmax)';r0=r;Nr=length(r);u=zeros(Nr,2);u(:,1)=u0;u(1,:)=ub;u(end,:)=ub;%%calcr=r(2:end-1);t=0;i=1;while(t=tmax)u(2:end-1,3-i)=u(2:end-1,i)+a*dt./(r.^2).*(r.^2.*(u(3:end,i)+u(1:end-2,i)-2*u(2:end-1,i))/dr^2+r.*(u(3:end,i)-u(1:end-2,i))./(2*dr)-u(2:end-1,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);endend参考书：顾昌鑫，计算物理学.上海：复旦大学出版社.2010