二维抛物型方程的隐式格式.docVIP

二维抛物型方程的隐式格式 2010-09-20 20:51 function varargout=liu(varargin)clear,clca=0;b=1;c=0;d=1;h1=1/30;h2=1/20;[X,Y,Z]=fc1t(a,b,c,d,h1,h2);figure()mesh(Y,X,Z);shading flat;xlabel(X,FontSize,14);ylabel(Y,FontSize,14);zlabel(error,FontSize,14);title(误差图)%求解下问题%u_xx+u_yy=u_t x,y 在区域内 axb,ayb%ctd%h1离散x方向的步长%h2离散t方向的步长function [X,Y,Z]=fc1t(a,b,c,d,h1,h2)x=a:h1:b;y=x;[X,Y]=meshgrid(x,y);t=c:h2:d;n=length(x);n1=length(t);%记录精确解UU=zeros(n,n,n1);for k=1:n1??? for r=1:n??????? for j=1:n??????????? UU(r,j,k)=fuu(x(r),y(j),t(k));??????? end??? endendN=100000;ee=10^(-8);s=h1^2/h2;%记录数值解 并赋初值U=zeros(n,n,n1);for i=1:n??? for j=1:n?????? U(i,j,1)=sef(x(i),y(j));??? endendfor r=1:n1??? for j=1:n?????? U(1,j,r)=seta3(y(j),t(r));??? endendfor r=1:n1??? for j=1:n?????? U(n,j,r)=seta4(y(j),t(r));??? endendfor r=1:n1??? for i=1:n?????? U(i,1,r)=seta1(x(i),t(r));??? endendfor r=1:n1??? for i=1:n?????? U(i,n,r)=seta2(x(i),t(r));??? endendfor k=2:n1??? for ki=1:N??????? U0=U(:,:,k);????? for i=2:n-1????????? for j=2:n-1??????????? U(i,j,k)=((U(i,j+1,k)+U(i,j-1,k)...??????????????? +U(i+1,j,k)+U(i-1,j,k))+s*U(i,j,k-1))/(s+4);???????? end????? end????? if max(max(abs(U0-U(:,:,k))))ee????????? break????? end??? endend Z=abs(UU(:,:,end)-U(:,:,end));function u=fuu(x,y,t)%精确解u=exp(0.5*(x+y)-t);function u=seta1(x,t)a=0;u=fuu(x,a,t);function u=seta2(x,t)b=1;u=fuu(x,b,t);function u=seta3(y,t)a=0;u=fuu(a,y,t);function u=seta4(y,t)%b=1b=1;u=fuu(b,y,t);function u=sef(x,y)%t=0u=fuu(x,y,0);