微分方程数值方法大作业.doc

微分方程数值方法大作业
一、 操作任务
给定初值问题
0<x<1,0<t<0.1
u(0,t)=u(1,t)=0 0<t<0.1
u(x,0)=sinπx 0≤x≤1
研究上述问题Grank-Nicholson格式的稳定性
二、 差分格式
先空间方向离散：
j=0,1,2……N

，，，，，，，
， ，，，，，，，，
取x=xj 得： ,
即得半离散格式：

j=0,1,2……N

记u(t)=[u(x1,t), u(x2,t) , …… u(xN-1,t)]T
F(t) =[f(x1,t)], f(x2,t)] …… f(xN-1,t)]T

格式变为：

U(0)=[]T
用梯形格式——Grank-Nicholason格式

则格式变为：

三、 程序
function x=EqtsForwardAndBackward(L,D,U,b)
%追赶法求解三对角线性方程组Ax=b
%x=EqtsForwardAndBackward(L,D,U,b)
%x:三对角线性方程组的解
%L:三对角矩阵的下对角线，行向量
%D:三对角矩阵的对角线，行向量
%U:三对角矩阵的上对角线，行向量
%b:线性方程组Ax=b中的b，列向量
n=length(D);m=length(b);
n1=length(L);n2=length(U);
if n-n1 ~= 1 || n-n2 ~= 1 || n ~= m %检查参数的输入是否正确
disp('输入参数有误！')
x=' ';
return;
end
%追的过程
for i=2:n
L(i-1)=L(i-1)/D(i-1);
D(i)=D(i)-L(i-1)*U(i-1);
end
x=zeros(n,1);
x(1)=b(1);
for i=2:n
x(i)=b(i)-L(i-1)*x(i-1);
end
%赶的过程
x(n)=x(n)/D(n);
for i=n-1:-1:1
x(i)=(x(i)-U(i)*x(i+1))/D(i);
end
return;
function [U x t] = PDEParabolicCN(uX,uT,M,N) %Crank-Nicolson隐式格式求解抛物型偏微分方程
%输入参数：uX -空间变量x的取值上限
% ? ?? ?? uT -时间变量t的取值上限
% ? ?? ?? M -沿x轴的等分区间数
% ? ?? ?? N -沿t轴的等分区间数
h=uX/M;%计算空间方向步长
tao=uT/N;%计算时间方向步长
x=(0:M)*h;
t=(0:N)*tao;
r=h/(tao*tao);%网格比
Diag=zeros(1,M-1);%矩阵的对角线元素
Low=zeros(1,M-2);%矩阵的下对角线元素
Up=zeros(1,M-2);%矩阵的上对角线元素
for i=1:M-2
Diag(i)=1+r;
Low(i)=-r/2;
Up(i)=-r/2;
end
Diag(M-1)=1+r;
%计算初值和边值
U=zeros(M+1,N+1);
for i=1:M+1
U(i,1)=sin(pi*x(i));
end
for j=1:N+1
U(1,j)=0;
U(M+1,j)=0;
end
B=zeros(M-1,M-1);
for i=1:M-2
B(i,i)=1-r;
B(i,i+1)=r/2;
B(i+1,i)=r/2;
end
B(M-1,M-1)=1-r;
%逐层求解，需要使用追赶法（调用函数EqtsForwardAndBackward）
for j=1:N
b1=zeros(M-1,1);
b1(1)=r*(U(1,j+1)+ U(1,j))/2;
b1(M-1)=r*(U(M+1,j+1)+U(M+1,j))/2;
b=B*U(2:M,j)+b1;
U(2:M,j+1)=EqtsForwardAndBackward(Low,Diag,Up,b);
end
U=U';
disp('请输入uX')
uX=input('uX=');
disp('请输入uT')
uT=input('uT=');
disp('请输入M')
M=input('M=');
disp('请输入N')
N=input('N=');
[U x t] = PDEParabolicCN(uX,uT,M,N);
mesh(x,t,U);
title('Grank-Nicholson，隐式格式，一维热传导方程解的图像')
xlabel('空间变量 x')
ylabel('时间变量t')
zlabel('一维热传导方程的解 x')
return;
四、 数据结果
M=50时：
M=100时：

五、 结论
该方程组Grank-Nicholason格式稳定.
六、 参考文献
1]　Garth A,Baker.Error estimates for finite element methods for second order hyperbolic equations[J].SIAM J Number Anal,1976,13(4):
564-576.
[2]　王申林.拟线性抛物型方程的Galerkin方法及H1模误差估计[J].山东大学学报,1989,24(1):19-31.
[3]　王自强,曹俊英,宋士仓.半线性抛物型方程改进全离散双尺度有限元分析[J].河南师范大学学报:自然科学版,2010,38(1):41-43.
[4]　Weiser A,Wheeler M F.On convergence of block-centered finite differences for elliptic problems[J].SIAM J Numer Anal,1982,25(2):
871-885.
[5]　王申林,孙淑英.对流—扩散问题的特征———块中心差分法[J].计算数学,1999,21(4):463-474.
5 / 5