最简差分格式

最简差分格式 朱勇军             20083738 考虑问题： 考虑一维热传导方程： ，            ， 其中a为正常数， 是给定的连续函数 边值条件： 网格剖分： 取空间步长 和时间步长 ，其中N，M都是正整数，用两族平行直线 和 将矩形域 分割成矩形网格，网格节点为 差分格式： 向前差分格式： 向量形式： 其中 ，以 表示网比，则该格式可变为： 矩阵形式： 向后差分格式： 向量形式： 其中 ，以 表示网比，则该格式可变为： 矩阵形式： -格式： 向量形式： 将向前差分格式和向后差分格式作加权平均，即得： 其中 ，以 表示网比，则该格式为： 矩阵形式： 特别当 时，即为下面的六点对称格式 六点对称格式： 向量形式： 将向前差分格式和向后差分格式作算术平均，即得： 其中 ，以 表示网比，则该格式为： 矩阵形式： 算例： 编程实现： 函数： 文件名：uxt.m function u=uxt(x,t) %求边值 if t==0 u=1-x; end if x==0 u=1; end if x==1 u=0; end %求真解 u=1-x; 文件名：fxt.m function f=fxt(x) % f=0; 向前差分格式程序： 文件名：forward.m clear; %输入区间信息,正常数a及时间和空间等分数 L=input('请输入空间上限L:'); T=input('请输入时间上限T:'); a=input('请输入正常数a:'); N=input('请输入空间等分数N:'); M=input('请输入时间等分数M:'); %计算时间和空间步长，网比r及节点值 h=L/N; s=T/M; r=a*s/(h*h); x=h:h:L; t=s:s:T; %合成系数矩阵A d1=r*ones(1,(N-2)); d2=(1-2*r)*ones(1,(N-1)); A=diag(d1,-1)+diag(d1,1)+diag(d2); %合成初值条件ux for p=1:N-1 ux(p)=uxt(x(p),0); end ux=ux'; %合成矩阵f for p=1:N-1 f(p)=s*fxt(x(p)); end f(1)=f(1)+r; f=f'; %计算u for q=1:M ux=A*ux+f; u(q,:)=ux'; ux=ux;    end disp('数值解：') u disp('精确解：') for q=1:M for p=1:N-1 ux(q,p)=uxt(x(p),0); end end ux 调试结果： >> forward 请输入空间上限L:1 请输入时间上限T:1 请输入正常数a:0.5 请输入空间等分数N:4 请输入时间等分数M:4 数值解： u = 0.7500    0.5000    0.2500 0.7500    0.5000    0.2500 0.7500    0.5000    0.2500 0.7500    0.5000    0.2500 精确解： ux = 0.7500    0.5000    0.2500 0.7500    0.5000    0.2500 0.7500    0.5000    0.2500 0.7500    0.5000    0.2500 >> 向后差分格式程序： 文件名：back.m %向后差分格式 clear; %输入区间信息,正常数a及时间和空间等分数 L=input('请输入空间上限L:'); T=input('请输入时间上限T:'); a=input('请输入正常数a:'); N=input('请输入空间等分数N:'); M=input('请输入时间等分数M:'); %计算时间和空间步长，网比r及节点值 h=L/N; s=T/M; r=a*s/(h*h); x=h:h:L; t=s:s:T; %合成系数矩阵A d1=-r*ones(1,(N-2)); d2=(1+2*r)*ones(1,(N-1)); A=diag(d1,-1)+diag(d1,1)+diag(d2); %合成初值条件ux for p=1:N-1 ux(p)=uxt(x(p),0); end ux=ux';    %合成矩阵f for p=1:N-1 f(p)=s*fxt(x(p)); end f(1)=f(1)+r; f=f'; %计算u for q=1:M ux=inv(A)*(ux+f); u(q,:)=ux'; ux=ux;    end disp('数值解：') u disp('精确解：') for q=1:M for p=1:N-1 ux(q,p)=uxt(x(p),0); end end ux 调试结果： >> back 请输入空间上限L:1 请输入时间上限T:1 请输入正常数a:0.5 请输入空间等分数N:5 请输入时间等分数M:5 数值解： u = 0.8000    0.6000    0.4000    0.2000 0.8000    0.6000    0.4000    0.2000 0.8000    0.6000    0.4000    0.2000 0.8000    0.6000    0.4000    0.2000 0.8000    0.6000    0.4000    0.2000 精确解： ux = 0.8000    0.6000    0.4000    0.2000 0.8000    0.6000    0.4000    0.2000 0.8000    0.6000    0.4000    0.2000 0.8000    0.6000    0.4000    0.2000 0.8000    0.6000    0.4000    0.2000 >> -格式程序： 文件名：seta.m %seta格式 clear; %输入区间信息,正常数a及时间和空间等分数 L=input('请输入空间上限L:'); T=input('请输入时间上限T:'); a=input('请输入正常数a:'); b=input('请输入权值b：'); N=input('请输入空间等分数N:'); M=input('请输入时间等分数M:'); %计算时间和空间步长，网比r及节点值 h=L/N; s=T/M; r=a*s/(h*h); x=h:h:L; t=s:s:T; %合成系数矩阵A d1=-(1-b)*r*ones(1,(N-2)); d2=(1+2*(1-b)*r)*ones(1,(N-1)); A=diag(d1,-1)+diag(d1,1)+diag(d2); %合成系数矩阵B b1=b*r*ones(1,(N-2)); b2=(1-2*b*r)*ones(1,(N-1)); B=diag(b1,-1)+diag(b1,1)+diag(b2); %合成初值条件ux for p=1:N-1 ux(p)=uxt(x(p),0); end ux=ux'; %合成矩阵f for p=1:N-1 f(p)=s*fxt(x(p)); end f(1)=f(1)+r; f=f'; %计算u for q=1:M ux=inv(A)*(B*ux+f); u(q,:)=ux'; ux=ux;    end disp('数值解：') u disp('精确解：') for q=1:M for p=1:N-1 ux(q,p)=uxt(x(p),0); end end Ux 调试结果： >> seta 请输入空间上限L:1 请输入时间上限T:1 请输入正常数a:0.5 请输入权值b：0.5 请输入空间等分数N:5 请输入时间等分数M:4 数值解： u = 0.8000    0.6000    0.4000    0.2000 0.8000    0.6000    0.4000    0.2000 0.8000    0.6000    0.4000    0.2000 0.8000    0.6000    0.4000    0.2000 精确解： ux = 0.8000    0.6000    0.4000    0.2000 0.8000    0.6000    0.4000    0.2000 0.8000    0.6000    0.4000    0.2000 0.8000    0.6000    0.4000    0.2000 >>