微分方程数值方法大作业1.docVIP

微分方程数值方法大作业 操作任务 给定初值问题 0x1,0t0.1 u(0,t)=u(1,t)=0 0t0.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(请输入u