偏微分方程的编程解法 2

[niedongbing]

前几天发了一个求助贴，不过没人回应，很失望，于是在失望中终于把这个方程搞定了，是高手们通常所说的挺简单的那种，不过没做出之前对我这个菜鸟来说挺难。现在把程序贴上供大家参考。（下面方程可调用MATLAB内部函数解决，但对建立起来的数学模型若不是标准的偏微分方程的话必须编程）
方程如下：dT2/dr2+(1/r)*dT/dr=dT/dt
I.C. : T=0
B.C. : r=1, dT/dr=-1/(2*pi)
r→∞, T=0
方程是一维的圆柱坐标方程
步骤如下：1.显示差分方程（r→∞, T=0 处理成r=1000,T=0） 2.编写程序
以下是程序：(建立两个M文件)
function [T,r,t] = NDB(D,tf,it0,brf,N,I)
% D为一向量，左右界限r0,rf是其元素
% tf为时间下限
% it0为函数T的初值
% brf分别为右边界值，左边界为牛顿边界用差分去得到
% N,I为r轴和时间轴区间数

r0 = D(1); rf = D(2);
dr = (rf-r0)/N; r = r0+[0:N]*dr;
dt = tf/I; t = [0:I]'*dt; %建立内节点


% I.C. 初值
for n = 1:N+1
T(n,1) = it0(r(n));
end
% B.C. 右边界值
for i = 1:I+1
T(N+1,i) = brf(t(i));
end

m1 = dt/dr/dr; m2 = dr/dt; m3 = dt/dr/pi;


% Solution 内节点的值（显式差分偏微分方程得到）
for i = 1:I
for n = 2:N
T(1,i+1) = (2*m1+1/m2)*T(2,i)+(1-2*m1-1/m2)*T(1,i)+ m3; % 左边界值，差分牛顿边值得到
T(n,i+1) = (m1+1/(m2+(n-1)/m1))*T(n+1,i)+(1-2*m1-1/(m2+(n-1)/m1))*T(n,i)+m1*T(n-1,i); %内节点值，显式差分方程得到

end
end



function zp()

r0 = 1; rf = 1000; %左边界为1，右边界为1000
N = 999; I = 4000; tf = 100; % r轴等分区间数为999，时间轴区间数为4000（可根据收敛性调整）
it0 = inline('0'); % 初始值

brf = inline('0'); %右边界值
D =[r0 rf];
[T,r,t] = NDB(D,tf,it0,brf,N,I);
% 以下是画图程序


figure
plot(r,T(:,40))
axis([0 100 0 0.1])
title('Solution')
xlabel('r')
ylabel('T')
hold on
plot(r,T(:,400))
hold on
plot(r,T(:,4000))
gtext('t=1'); gtext('t=10'); gtext('t=100')

上面由于采用显式差分，存在收敛性问题，上面并未写入判断收敛与否的程序，而是通过合理设置N和I的值来保证收敛的。调用内部函数的解法一般参考书上都有的，就不贴了。