m=1/hi;n=t/ti;
g=ti/(hi^2); %--------------------------------- g为网格比
u=ones(m+1,m+1); %------------------------输入初始值
for i=2:m
for j=2:m
u(i,j)=sin(pi*(i-1)*hi)*sin(2*pi*(j-1)*hi);
1.00000000000000
1.00000000000000
0.99999999970799
0.99999999944526
0.99999999923550
0.99999999910241
0.99999999905502
0.99999999910240
0.99999999923550
0.99999999944526
0.99999999970799
1.00000000000000
1.00000000000000
0.99999999944526
0.99999999894348
0.99999999854766
0.99999999829052
0.99999999820481
0.99999999829052
0.99999999854766
function pr(ti,h,t)
%-------------------------------------------ti:时间步长h:空间步长;
k=t/ti+1;
m=1/h+1;
r=ti/h^2; %------------------------------ r为网格比;
w=ones(m,m);
end
end
a(1:m-2)=-0.5*g;
b(1:m-1)=1+g;
c(1:m-2)=-0.5*g; %------------------------输入用局部一维差分格式求解的三对角矩阵
B=zeros(m-1,m+1);
for i=1:m-1
B(i,i)=0.5*g;
B(i,i+1)=1-g;
u(i,j)=sin(pi*(i-1)*h)*sin(2*pi*h*(j-1));
end
end
tic
for l=1:k
for i=2:m-1
for j=2:m-1
w(i,j)=r*u(i-1,j)+r*u(i,j-1)+r*u(i+1,j)+r*u(i,j+1)+(1-4*r)*u(i,j);
end
u=ones(m,m); %------------------------输入初始值
v=ones(m,m);
for i=2:m-1
for j=2:m-1
u(i,j)=sin(pi*(i-1)*h)*sin(2*pi*h*(j-1));
end
end
%------------------------输入用P-R差分格式求解的三对角矩阵
0.99999999694199
0.99999999709532
0.99999999752602
0.99999999820481
0.99999999905502
1.00000000000000
1.00000000000000
0.99999999910240
0.99999999829052
0.99999999765007 >'>
p(1)=2*r;
p(m-2)=2*r;
tic
for l=1:k
for i=2:m-1
d1=A*u(i,2:m-1)'+p;
d1=d1';
w(2:m-1,i)=zg(a,b,c,d1); %-------------------------调用追赶法求解
d2=A*w(2:m-1,i)+p;
v(i,2:m-1)=zg(a,b,c,d2); %-------------------------调用追赶法求解
温馨提示：如果这篇偏微分方程数值解上机文档不完整，请从下面的
0.99999999799851
0.99999999854766
0.99999999923550
1.00000000000000
1.00000000000000
0.9999999991ቤተ መጻሕፍቲ ባይዱ240
0.99999999829052
0.99999999765007
0.99999999723401
0.99999999709532
end
x(n)=y(n)/u(n); %------------------------追赶法求解之赶过程求解Uz=y;
for j=n-1:-1:1
if u(j)==0
break;
else
x(j)=(y(j)-c(j)*x(j+1))/u(j);
end
end
%-----------------------------------------------运用P-R差分格式求解二维扩散方程的初边值问题;
%--------------------------------------------d:追赶法所求方程的右端向量;
%--------------------------------------------l:系数矩阵A所分解成的下三角阵L中的下对角元素了l(i);
%--------------------------------------------u:系数矩阵A所分解成的下三角阵U中的主对角元素了u(i);
%--------------------------------------------a:方程组系数矩阵A的下对角元素;
%--------------------------------------------b:方程组系数矩阵A的主对角元素;
%--------------------------------------------c:方程组系数矩阵A的上对角元素;
end
u=v';
end
toc
t=toc
u
mesh(0:0.1:1,0:0.1:1,u)
局部一维格式:
将原格式化为:
代入边界条件,转化为三对角矩阵
附源程序:
%------------------------------------------运用局部一维格式求解二维扩散方程的初边值问题;
function god(ti,hi,t) %------------------------------------------------ti为时间步长, hi为空间步长;
1.00000000000000
1.00000000000000
1.00000000000000
1.00000000000000
1.00000000000000
1.00000000000000
1.00000000000000
1.00000000000000
1.00000000000000
1.00000000000000
end
u=w;
end
toc
t=toc
u
mesh(u)
交替方向隐式格式(P-R格式):
将原差分格式化为:
代入边界条件,转化为三对角矩阵
附追赶法源程序:
%-------------------------------------------追赶法求解三对角方程组;
function x=zg(a,b,c,d)
n=length(b);
u(1)=b(1);
y(1)=d(1);
for i=1:n-1 %--------------------------追赶法求解之追过程求解Ly=d;
l(i)=a(i)/u(i);
u(i+1)=b(i+1)-l(i)*c(i);
y(i+1)=d(i+1)-l(i)*y(i);
0.99999999723401
0.99999999765007
0.99999999829052
0.99999999910240
1.00000000000000
1.00000000000000
0.99999999905502
0.99999999820481
0.99999999752602
0.99999999709532
B(i,i+2)=0.5*g;
end
f=zeros(m-1,1);
f(1,1)=0.5*g;
f(m-1,1)=0.5*g;
w=ones(m+1,m+1);
for i=1:n
for j=2:m
d=B*u(:,j)+f;
%-------------------------调用追赶法求解
x=zg(a,b,c,d);
b=ones(1,m-2)*(2+2*r);
a=-r*ones(1,m-3);
c=-r*ones(1,m-3);
A=zeros(m-2,m-2);
for i=1:m-2
A(i,i)=2-2*r;
end
for i=1:m-3
A(i,i+1)=r;
A(i+1,i)=r;
end
p=zeros(m-2,1);
数值求解二维扩散方程的初边值问题
古典显式格式:
将原格式化为:
附源程序:
%-------------------------------------------运用古典显式差分格式求解二维扩散方程的初边值问题;