ADI(交替隐式求求解抛物型偏微分方程)程序.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % Example of ADI Method for 2D heat equation % % % % u_t = u_{xx} + u_{yy} + f(x,t) % %a<x<b;c<y<d % % Test problme: % % Exact solution: u(t,x,y) = exp(-t) sin(pi*x)sin(pi*y) % % Source term: f(t,x,y) = exp(-t) sin(pi*x) sin(pi*y) (2pi^2-1) % % % % Files needed for the test: % % % % adi.m: This file, the main calling code. % % f.m: The file defines the f(t,x,y) % % uexact.m: The exactsolution. % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear; close all; a = 0; b=1; c=0; d=1; n = 40; tfinal = 1; m = n; h = (b-a)/n; dt=h; h1 = h*h; r=dt/h1; x=a:h:b; y=c:h:d; %-- Initial condition: t = 0; for i=1:m+1, for j=1:m+1, u1(i,j) = uexact(t,x(i),y(j)); end end %---------- Big loop for time t -------------------------------------- k_t = fix(tfinal/dt); for k=1:k_t t1 = t + dt; t2 = t + dt/2; %--- sweep in x-direction -------------------------------------- for i=2:m, % Boundary condition. u2(1,i) = -r/2*uexact(t1,x(1),y(i-1)) + (1+r)*uexact(t1,x(1),y(i)) + -r/2*uexact(t1,x(1),y(i+1));u2(m+1,i) = -r/2*uexact(t1,x(m+1),y(i-1)) + (1+r)*uexact(t1,x(m+1),y(i)) + -r/2*uexact(t1,x(m+1),y(i+1)); end for j = 2:n, % Look for fixed y(j) b=zeros(m-1,1); for i=2:m, b(i-1) = r^2/4*( u1(i-1,j-1) + u1(i-1,j+1) + u1(i+1,j-1) + u1(i+1,j+1) ) ... +r/2*( 1-r )*( u1(i-1,j) + u1(i+1,j) + u1(i,j-1) + u1(i,j+1) )... + ( 1-2*r+r^2 )*u1(i,j) -3/2*dt*exp( 1/2*( (i-1)*h + (j-1)*h ) -t2 ); end b(1) = b(1) + r/2 * u2(1,j); b(m-1) =b(m-1) + r/2 * u2(m+1,j); A = diag( (1+r)*ones(m-1,1) ) + diag( -r/2*ones(m-2,1),1 ) ... + diag(-r/2*ones(m-2,1),-1); ut = A\b; % Solve the diagonal matrix. for i=1:m-1,u2(i+1,j) = ut(i); end end % Finish x-sweep. %-------------- loop in y -direction -------------------------------- for i=1:m+1, % Boundary condition u1(i,1) = uexact(t1,x(i),y(1));u1(i,n+1) = uexact(t1,x(i),y(m+1)); u1(1,i) = uexact(t1,x(1),y(i)); u1(m+1,i) =uexact(t1,x(m+1),y(i)); end for i = 2:m, for j=2:m, b(j-1) = u2(i,j); end b(1) = b(1) + r/2*u1(i,1); b(m-1) = b(m-1) + r/2*u1(i,m+1);A = diag( (1+r)*ones(m-1,1) ) + diag( -r/2*ones(m-2,1),1 ) ... + diag(-r/2*ones(m-2,1),-1); ut = A\b; for j=1:n-1, u1(i,j+1) = ut(j); end end % Finish y-sweep. t = t + dt; %--- finish ADI method at this time level, go to the next time level. end %-- Finished with the loop in time %----------- Data analysis ---------------------------------- for i=1:m+1, for j=1:n+1, ue(i,j) = uexact(tfinal,x(i),y(j)); end end [xx,yy]=meshgrid(x,y); xx=xx';yy=yy'; mesh(xx,yy,u1) title('数值解图') pause surf(xx,yy,ue) title('精确解图') pausesurf(xx,yy,u1-ue) title('误差图') e = max(max(abs(u1-ue))) % The infinityerror. %******%表示原来程序编写错误 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % Example of ADI Method for 2D heat equation % % % % u_tt = u_{xx} + u_{yy} + f(x,t) % %a<x<b;c<y<d % % Test problme: % % 精确解: u(t,x,y) = exp( 1/2*(x+y) - t ) % % 右端函数: f(t,x,y) = 1/2*exp( 1/2*(x+y) - t ) % % 初始条件: g1 = exp( 1/2*(x+y) ) % % 初始导数条件: g2 = -exp( 1/2*(x+y) ) % % 需要函数 % % g1.m g2.m 初始条件 % % f2.m: The file defines the f(t,x,y) % % uexact.m: The exactsolution. % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear; close all; a = 0; b=1; c=0; d=1; n = 10; tfinal = 1; m = n; h = (b-a)/n; dt=h; s=(dt/h)^2; x=a:h:b; y=c:h:d; %-- Initial condition: t = 0; for i=1:m+1, for j=1:m+1, u1(i,j) = g1(x(i),y(j)); end end %-- Initial deriavtive condition: for i=1:m+1, for j=1:m+1, %******% u2(i,j) = g1(x(i),y(j)) + dt * g2( x(i),y(j) ) + dt^2/4 * ( g1(x(i),y(j)) + f2 (t, x(i),y(j) ) ); u2(i,j) =g1(x(i),y(j)) + dt * g2( x(i),y(j) ) + dt^2/2 * ( g1(x(i),y(j))/2 + f2 (t, x(i),y(j) ) ); endend %---------- Big loop for time t -------------------------------------- k_t = fix(tfinal/dt); for k=2:k_t t1 = t + dt; t2 = t1 + dt; %--- sweep in x-direction --------------------------------------for i=1:m+1, % Boundary condition. u4(1,i) = ( uexact(t,x(1),y(i)) -2*uexact(t1,x(1),y(i)) + uexact(t2,x(1),y(i)) )/dt^2 ; u4(m+1,i) = ( uexact(t,x(m+1),y(i)) - 2*uexact(t1,x(m+1),y(i)) + uexact(t2,x(m+1),y(i)) )/dt^2 ; u4(i,1) = ( uexact(t,x(i),y(1)) - 2*uexact(t1,x(i),y(1)) + uexact(t2,x(i),y(1)) )/dt^2 ; u4(i,m+1) = ( uexact(t,x(i),y(m+1)) - 2*uexact(t1,x(i),y(m+1)) + uexact(t2,x(i),y(m+1)) )/dt^2 ; end for i=2:m u3(1,i)=-s/2*u4(1,i-1) + (1+s)*u4(1,i) - s/2*u4(1,i+1); u3(m+1,i)=-s/2*u4(m+1,i-1) +(1+s)*u4(m+1,i) - s/2*u4(m+1,i+1); end for j = 2:n, % Look for fixed y(j) b=zeros(m-1,1); for i=2:m, b(i-1) = ( u2(i-1,j) + u2(i+1,j) + u2(i,j-1) + u2(i,j+1) -4*u2(i,j) )/h^2 +f2 (t1, x(i),y(j) ); end b(1) = b(1) + s/2 * u3(1,j); b(m-1) = b(m-1) + s/2 * u3(m+1,j); A = diag( (1+s)*ones(m-1,1) ) + diag( -s/2*ones(m-2,1),1 ) ... + diag(-s/2*ones(m-2,1),-1); ut = A\b; % Solve the diagonal matrix. for i=1:m-1, u3(i+1,j) = ut(i); end end % Finish x-sweep. %-------------- loop in y -direction -------------------------------- for i = 2:m, forj=2:m, b(j-1) = u3(i,j); end b(1) = b(1) + s/2*u4(i,1); b(m-1) = b(m-1) + s/2*u4(i,m+1); ut = A\b; for j=1:n-1, u4(i,j+1) = ut(j); end end % Finish y-sweep. for i=1:m+1, % Boundary condition u5(i,1) = uexact(t2,x(i),y(1)); u5(i,m+1) = uexact(t2,x(i),y(m+1));u5(1,i) = uexact(t2,x(1),y(i)); u5(m+1,i) = uexact(t2,x(m+1),y(i)); end for i=2:m, forj=2:m, u5(i,j) = dt^2 * u4(i,j) + 2*u2(i,j) - u1(i,j); end end u1=u2; u2=u5; %******%t = t2 + dt; t = t+dt ; %--- finish ADI method at this time level, go to the next time level. end %-- Finished with the loop in time %----------- Data analysis ---------------------------------- for i=1:m+1, for j=1:n+1, ue(i,j) = uexact(tfinal,x(i),y(j)); end end[xx,yy]=meshgrid(x,y); xx=xx';yy=yy'; mesh(xx,yy,u2) title('数值解图') pausesurf(xx,yy,ue) title('精确解图') pause %******%surf(xx,yy,u2-ue) surf(xx,yy,u2-ue) title('误差图') e = max(max(abs(u2-ue))) %******%e = max(max(abs(u1-ue))) % The infinityerror. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % Example of ADI Method for 2D heat equation % % % % u_t = u_{xx} + u_{yy} + f(x,t) % %a<x<b;c<y<d % % Test problme: % % Exact solution: u(t,x,y) = exp(-t) sin(pi*x) sin(pi*y) % % Source term: f(t,x,y) = exp(-t) sin(pi*x) sin(pi*y) (2pi^2-1) % % % % Files needed for the test: % % % % f.m: The file defines the f(t,x,y) % % uexact.m: The exact solution. % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear; close all; a = 0; b=1; c=0; d=1; n = 20; tfinal = 1; m = n; h = (b-a)/n; h1 = h*h; dt=h1; r=dt/h1; x=a:h:b; y=c:h:d; %-- Initial condition: t = 0; for i=1:m+1, for j=1:m+1, u1(i,j) = uexact(t,x(i),y(j)); end end %---------- Big loop for time t -------------------------------------- k_t = ceil(tfinal/dt); for k=1:k_t t1 = t + dt; t2 = t + dt/2; %--- sweep in x-direction -------------------------------------- for j=2:m, % Boundary condition. u2(1,j) = (1/12-r/2)*uexact(t1,x(1),y(j-1)) + (5/6+r)*uexact(t1,x(1),y(j)) + (1/12-r/2) *uexact(t1,x(1),y(j+1)); u2(m+1,j) = (1/12-r/2)*uexact(t1,x(m+1),y(j-1)) +(5/6+r)*uexact(t1,x(m+1),y(j)) + (1/12-r/2) * uexact(t1,x(m+1),y(j+1)); end for j =2:n, % Look for fixed y(j) b=zeros(m-1,1); for i=2:m, b(i-1) = (1/144+r/12+r^2/4) * ( u1(i-1,j-1) + u1(i-1,j+1) + u1(i+1,j-1) + u1(i+1,j+1) ) ... + ( 10/144+r/3-r^2/2 )*( u1(i-1,j) + u1(i+1,j) + u1(i,j-1) + u1(i,j+1) )... + ( 100/144 - 40 * r/24 + r^2 )*u1(i,j)... + dt * ( f(t2,x(i-1),y(j-1)) + f(t2,x(i-1),y(j+1)) + f(t2,x(i+1),y(j-1)) + f(t2,x(i+1),y(j+1))... + 10 * ( f(t2,x(i-1),y(j)) + f(t2,x(i+1),y(j)) + f(t2,x(i),y(j-1)) + f(t2,x(i),y(j+1)) )... + 100 *f(t2,x(i),y(j)) ) / 144; end b(1) = b(1) + ( r/2 -1/12 ) * u2(1,j); b(m-1) = b(m-1) + ( r/2 -1/12 ) * u2(m+1,j); A = diag( (5/6+r)*ones(m-1,1) ) + diag( (1/12-r/2) * ones(m-2,1),1 ) ... + diag( (1/12-r/2) * ones(m-2,1),-1); ut = A\b; % Solve the diagonal matrix. for i=1:m-1, u2(i+1,j) = ut(i); end end % Finish x-sweep. %-------------- loop in y -direction -------------------------------- for i=1:m+1, % Boundary condition u1(i,1) = uexact(t1,x(i),y(1)); u1(i,n+1) = uexact(t1,x(i),y(m+1)); u1(1,i) = uexact(t1,x(1),y(i));u1(m+1,i) = uexact(t1,x(m+1),y(i)); end for i = 2:m, for j=2:m, b(j-1) = u2(i,j); end b(1) =b(1) + ( r/2 -1/12 ) * u1(i,1); b(m-1) = b(m-1) + ( r/2 -1/12 ) * u1(i,m+1); ut = A\b; for j=1:n-1, u1(i,j+1) = ut(j); end end % Finish y-sweep. t = t + dt; %--- finish ADI method at this time level, go to the next time level. end %-- Finished with the loop in time %----------- Data analysis ---------------------------------- for i=1:m+1, for j=1:n+1, ue(i,j) = uexact(tfinal,x(i),y(j)); end end [xx,yy]=meshgrid(x,y); xx=xx';yy=yy'; mesh(xx,yy,u1) title('数值解图') pause surf(xx,yy,ue) title('精确解图') pause surf(xx,yy,u1-ue) title('误差图') e = max(max(abs(u1-ue))) % The infinity error. %算例1.1.3，见课本第17页 %首先网格剖分 xl=0.; xr=pi;%左右端点 n=20; h=(xr-xl)/n;%步长 %解差分格式形成的方程组 U=zeros(n-1,1); aa=-1+h^2/12; bb=2+5/6*h^2; A=diag(bb*ones(n-1,1))+diag(aa*ones(n-2,1),1)+diag(aa*ones(n-2,1),-1);%系数矩阵 F=zeros(n-1,1); for i=1:n-1 F(i)=h^2/12*(exp(xl+(i-1)*h)*(sin(xl+(i-1)*h)-2*cos(xl+(i-1)*h))+10*exp(xl+i*h)*(sin(xl+i*h)-2*cos(xl+i*h))+exp(xl+(i+1)*h)*(sin(xl+(i+1)*h)-2*cos(xl+(i+1)*h)));%右端项 end U=inv(A)*F;%得到数值解 %精确解与数值解误差绝对值的最大值 ana=zeros(n-1,1); for i=1:n-1 ana(i)=exp(xl+i*h)*sin(xl+i*h);%精确解 end wucha=max(abs(U-ana)) %可视化处理 x=xl:h:xr; U=[0 U' 0];%U表示数值解ana=[0 ana' 0];%ana表示精确解 plot(x,U,'r*',x,ana); %向后Euler差分格式 %Ut-Uxx=0,0<x<1,0<t<1 %u(x,0)=exp(x) %u(0,t)=exp(t),u(1,t)=exp(1+t) %%%%%%%网格剖分%%%%%%%%% % first set up the attributes h = 0.1; % the step-size in the x dimension%空间步长 r=1; %步长比 dt = h^2 * r; % the step-size in the t dimension%%时间步长 T0 = 1; % the time-interval over which we will calculate u 时间长度 TN = ceil(T0 / dt); % the number of steps in the time-interval%时间网格 X0 = 1; % the x-interval over which we will calculate u空间长度 XN = ceil(X0 / h); % the number of steps in the x-interval%空间网格 U = zeros(XN + 1, TN + 1); % set up an array of zeros %初始化 % Next we set up the initial conditions%初始值及边界值 x = 0 :h: X0; for j=1:XN+1 U(j , 1) = exp( (j-1)*h ); end for k=1:TN+1 U(1,k) =exp( (k-1) *dt);U(XN+1,k)=exp(1+ ( k-1 ) *dt); end %%%%%%形成系数矩阵%%%%%%%% A=diag((5/6+r)*ones(XN-1,1))+diag( (1/12-r/2) *ones(XN-2,1),1 )+diag( (1/12-r/2)*ones(XN-2,1),-1 ); B=diag((5/6-r)*ones(XN-1,1))+diag( (1/12+r/2) *ones(XN-2,1),1 )+diag( (1/12+r/2)*ones(XN-2,1),-1 ); %右端项初始化 F = zeros( XN-1,1 ); %迭代求解 for k = 1 : TN F(1)=(1/12+r/2)*exp((k-1)*dt)-(1/12-r/2)*exp(k*dt); F(XN-1)=(1/12+r/2)*exp((k-1)*dt+1)-(1/12-r/2)*exp(k*dt+1); U(2:XN,k+1)= A \ ( B * U(2:XN,k) + F ); end; %给出精确解 ana = zeros(XN + 1, TN + 1); for k = 1 : TN+1 for j = 1 : XN+1 ana(j, k ) = exp ( (j-1) *h + (k-1)*dt); end; end; %给出在T=1时的最大误差 max( abs ( U (:,TN+1) -ana (:,TN+1) ) ) %图示化 t=0:dt:T0; [xx,tt]= meshgrid(x,t); xx=xx'; tt=tt'; mesh(xx,tt,ana); title('analytical solution');ylabel(sprintf('steps in x domain\nfrom 0 to 1')); xlabel(sprintf('steps in time-domain\nfrom 0 to 1')); zlabel('U(x,y)'); pause mesh(xx,tt,U); title('numerical solution'); ylabel(sprintf('steps in x domain\nfrom 0 to 1')); xlabel(sprintf('steps in time-domain\nfrom 0 to 1')); zlabel('U(xi,tk)'); pause mesh(xx,tt,U-ana); title('diffeence beween analical and numeical soluion'); ylabel(sprintf('steps in x domain\nfrom 0 to 1')); xlabel(sprintf('steps in time-domain\nfrom 0 to 1')); zlabel('U(xi,tk)'); pauseplot(x,U(:,TN+1),'-or',x,ana(:,TN+1)) title('向前Euler格式数值解与精确解的比较') legend('uh(x,1)','u(x,1)')%Crank-Nicolson schemes % 显差分格式 %算例4.1.1 %%%%%%%网格剖分%%%%%%%%% h = 1/10; %%空间步长 a = 1; dt =1/20; %时间步长 s = a * dt/h; %步长比 T = 1; % 时间长度 n = ceil(T / dt); % 时间网格 X0 = 1; % 空间长度 m = ceil(X0 / h); %空间网格 u = zeros(n + 1, m + 1); %%数解初始化u(tk,xi) % %初始值及边界值 x=0:h:X0; phi=exp(x); u(1,:)=phi; %u的第一行表示初始状态 t=0:dt:T; u(:,1)=exp(t'); %左界u(:,m+1)=exp(1+t'); %%%%%%%%%%%%%%% f=zeros(n-1,m-1); for k=1:n+1 for i=1:m+1 f(k,i)=0; end end %%%%%%%%%%%%%%% psi=exp(x); for i = 2:m u(2,i) = dt * psi(i) + s^2/2 * ( phi(i-1) + phi(i+1) ) ... + ( 1 - s^2 ) * phi(i) + dt^2/2 * f(1,i);end %%%%%%%%%%%%%%% for k=2:n for i=2:m u (k+1 , i) = s^2 * u (k , i-1) +2*( 1-s^2 )*u (k , i) ... + s^2 * u (k , i+1) - u (k-1 , i) + dt^2 * f( k , i); endend %%%%%%%%%%%%%%%%给出精确解 ana = zeros(n + 1, m + 1); for k = 1 : n+1 for i = 1 : m+1 ana( k, i ) = exp ( (i-1) *h + (k-1)*dt); end;end; %%%%%%%%%%%%%%%给出在T=1时的最大误差 max( abs ( u (n+1,:) -ana (n+1,:) ) ) %图示化 [tt,xx]= meshgrid(t,x); xx=xx'; tt=tt'; mesh(tt,xx,ana);title('analytical solution'); xlabel(sprintf('steps in time domain\nfrom 0 to 1'));ylabel(sprintf('steps in x-domain\nfrom 0 to 1')); zlabel('U(x,y)'); pause mesh(tt,xx,u); title('numerical solution'); xlabel(sprintf('steps in time domain\nfrom 0 to 1'));ylabel(sprintf('steps in x-domain\nfrom 0 to 1')); zlabel('u(xi,tk)'); pause mesh(tt,xx,u-ana); title('diffeence beween analical and numeical soluion'); xlabel(sprintf('steps in time- domain\nfrom 0 to 1')); ylabel(sprintf('steps in x domain\nfrom 0 to 1')); zlabel('u-ana'); pause plot(x,u(n+1,:),'-or',x,ana(n+1,:)) title('显格式数值解与精确解的比较') legend('uh(x,1)','u(x,1)') % compace 差分格式 %算例4.3.1 %%%%%%%网格剖分%%%%%%%%% h = 1/10; %%空间步长 a = 1; dt =1/100; %时间步长 s = a * dt/h; %步长比 T = 1; % 时间长度 n = ceil(T / dt); % 时间网格 X0 = 1; % 空间长度 m = ceil(X0 / h); %空间网格 u = zeros(n + 1, m + 1); %%数解初始化u(tk,xi) % %初始值及边界值 x=0:h:X0; phi=exp(x); u(1,:)=phi; %u的第一行表示初始状态 t=0:dt:T; u(:,1)=exp(t'); %左界u(:,m+1)=exp(1+t'); %%%%%%%%%%%%%%% f=zeros(n+1,m+1); for k=1:n+1 fori=1:m+1 f(k,i)=0; end end %%%%%%%%%%%%%%% psi=exp(x); for i = 2:m u(2,i) = dt * psi(i) + s^2/2 * ( phi(i-1) + phi(i+1) ) ... + ( 1 - s^2 ) * phi(i) + dt^2/2 * f(1,i);end %%%%%%%%%%%%%%% A=diag( (5/6+s^2 )*ones( m-1,1 ) )+diag( (1/12-1/2*s^2) * ones (m-2,1),1 )+diag( (1/12-1/2*s^2) * ones (m-2,1),-1 ); B=diag( -(5/6+s^2 )*ones( m-1,1 ) )+diag( (1/2*s^2-1/12) * ones (m-2,1),1 )+diag( (1/2*s^2-1/12) * ones (m-2,1),-1 ); C=diag(5/3*ones(m-1,1))+diag(1/6*ones(m-2,1),1)+diag(1/6*ones(m-2,1),-1); for k=2:n for i=2:m F(i-1)=1/12*dt^2*( f(i-1)+10*f(i)+f(i+1) ); end F(1)=F(1) + (1/2*s^2-1/12) * ( u(k+1,1) + u(k-1,1) )+1/6*u(k,1); F(m-1)=F(m-1) + (1/2*s^2-1/12) * ( u(k+1,m+1) + u(k-1,m+1) )+1/6*u(k,m+1); y = A \ ( B * u(k-1,2:m )' + F' + C* u( k,2:m )' ); for i=1:m-1u(k+1,i+1)=y(i); end end %%%%%%%%%%%%%%%%给出精确解 ana = zeros(n + 1, m + 1); for k = 1 : n+1 for i = 1 : m+1 ana( k, i ) = exp ( (i-1) *h + (k-1)*dt); end; end; %%%%%%%%%%%%%%%给出在T=1时的最大误差 max( abs ( u (n+1,:) -ana (n+1,:) ) ) %图示化 [tt,xx]= meshgrid(t,x); xx=xx'; tt=tt'; mesh(tt,xx,ana);title('analytical solution'); xlabel(sprintf('steps in time domain\nfrom 0 to 1'));ylabel(sprintf('steps in x-domain\nfrom 0 to 1')); zlabel('U(x,y)'); pause mesh(tt,xx,u); title('numerical solution'); xlabel(sprintf('steps in time domain\nfrom 0 to 1'));ylabel(sprintf('steps in x-domain\nfrom 0 to 1')); zlabel('u(xi,tk)'); pause mesh(tt,xx,u-ana); title('diffeence beween analical and numeical soluion'); xlabel(sprintf('steps in time- domain\nfrom 0 to 1')); ylabel(sprintf('steps in x domain\nfrom 0 to 1')); zlabel('u-ana'); pause plot(x,u(n+1,:),'-or',x,ana(n+1,:)) title('显格式数值解与精确解的比较') legend('uh(x,1)','u(x,1)') % implicit difference scheme %算例4.2.1 %%%%%%%网格剖分%%%%%%%%% h = 1/10; %%空间步长 a = 1; dt = 1/10; %时间步长 s = a *dt/h; %步长比 T = 1; % 时间长度 n = ceil(T / dt); % 时间网格 X0 = 1; % 空间长度 m = ceil(X0 / h); %空间网格 u = zeros(n + 1, m + 1); %%数解初始化u(tk,xi) % %初始值及边界值 x=0:h:X0; phi=exp(x); u(1,:)=phi; %u的第一行表示初始状态 t=0:dt:T; u(:,1)=exp(t'); %左界u(:,m+1)=exp(1+t'); %%%%%%%%%%%%%%% f=zeros(n+1,m+1); for k=1:n+1 for i=1:m+1 f(k,i)=0; end end %%%%%%%%%%%%%%% psi=exp(x); for i = 2:m u(2,i) = dt * psi(i) + s^2/2 * ( phi(i-1) + phi(i+1) ) ... + ( 1 - s^2 ) * phi(i) + dt^2/2 * f(1,i);end %%%%%%%%%%%%%%% A=diag( (1+s^2 )*ones( m-1,1 ) )+diag( -1/2*s^2 * ones (m-2,1),1 )+diag( -1/2*s^2 * ones (m-2,1),-1 ); B=diag( -(1+s^2 )*ones( m-1,1 ) )+diag( 1/2*s^2 * ones (m-2,1),1 )+diag( 1/2*s^2 * ones (m-2,1),-1 ); for k=2:n F = dt^2 * f( k,2:m ); F(1)=F(1) + 1/2*s^2*( u(k+1,1) + u(k-1,1) ); F(m-1)=F(m-1) +1/2*s^2*( u(k+1,m+1) + u(k-1,m+1) ); y = A \ ( B * u(k-1,2:m )' + F' + 2 * u( k,2:m )' ); for i=1:m-1 u(k+1,i+1)=y(i); end end %%%%%%%%%%%%%%%%给出精确解 ana= zeros(n + 1, m + 1); for k = 1 : n+1 for i = 1 : m+1 ana( k, i ) = exp ( (i-1) *h + (k-1)*dt); end; end; %%%%%%%%%%%%%%%给出在T=1时的最大误差 max( abs ( u (n+1,:) -ana (n+1,:) ) ) %图示化 [tt,xx]= meshgrid(t,x); xx=xx'; tt=tt';mesh(tt,xx,ana); title('analytical solution'); xlabel(sprintf('steps in time domain\nfrom 0 to 1')); ylabel(sprintf('steps in x-domain\nfrom 0 to 1')); zlabel('U(x,y)'); pausemesh(tt,xx,u); title('numerical solution'); xlabel(sprintf('steps in time domain\nfrom 0 to 1')); ylabel(sprintf('steps in x-domain\nfrom 0 to 1')); zlabel('u(xi,tk)'); pausemesh(tt,xx,u-ana); title('diffeence beween analical and numeical soluion');xlabel(sprintf('steps in time- domain\nfrom 0 to 1')); ylabel(sprintf('steps in xdomain\nfrom 0 to 1')); zlabel('u-ana'); pause plot(x,u(n+1,:),'-or',x,ana(n+1,:)) title('显格式数值解与精确解的比较') legend('uh(x,1)','u(x,1)') % compace 差分格式 %算例4.3.1 %%%%%%%网格剖分%%%%%%%%% h = 1/10; %%空间步长 a = 1; dt =1/100; %时间步长 s = a * dt/h; %步长比 T = 1; % 时间长度 n = ceil(T / dt); % 时间网格 X0 = 1; % 空间长度 m = ceil(X0 / h); %空间网格 u = zeros(n + 1, m + 1); %%数解初始化u(tk,xi) % %初始值及边界值 x=0:h:X0; phi=exp(x); u(1,:)=phi; %u的第一行表示初始状态 t=0:dt:T; u(:,1)=exp(t'); %左界u(:,m+1)=exp(1+t'); %%%%%%%%%%%%%%% f=zeros(n+1,m+1); for k=1:n+1 for i=1:m+1 f(k,i)=0; end end %%%%%%%%%%%%%%% psi=exp(x); for i = 2:m u(2,i) = dt * psi(i) + s^2/2 * ( phi(i-1) + phi(i+1) ) ... + ( 1 - s^2 ) * phi(i) + dt^2/2 * f(1,i);end %%%%%%%%%%%%%%% A=diag( (5/6+s^2 )*ones( m-1,1 ) )+diag( (1/12-1/2*s^2) * ones (m-2,1),1 )+diag( (1/12-1/2*s^2) * ones (m-2,1),-1 ); B=diag( -(5/6+s^2 )*ones( m-1,1 ) )+diag( (1/2*s^2-1/12) * ones (m-2,1),1 )+diag( (1/2*s^2-1/12) * ones (m-2,1),-1 ); C=diag(5/3*ones(m-1,1))+diag(1/6*ones(m-2,1),1)+diag(1/6*ones(m-2,1),-1); for k=2:n for i=2:m F(i-1)=1/12*dt^2*( f(i-1)+10*f(i)+f(i+1) ); end F(1)=F(1) + (1/2*s^2-1/12) * ( u(k+1,1) + u(k-1,1) )+1/6*u(k,1); F(m-1)=F(m-1) + (1/2*s^2-1/12) * ( u(k+1,m+1) + u(k-1,m+1) )+1/6*u(k,m+1); y = A \ ( B * u(k-1,2:m )' + F' + C* u( k,2:m )' ); for i=1:m-1u(k+1,i+1)=y(i); end end %%%%%%%%%%%%%%%%给出精确解 ana = zeros(n + 1, m + 1); for k = 1 : n+1 for i = 1 : m+1 ana( k, i ) = exp ( (i-1) *h + (k-1)*dt); end; end; %%%%%%%%%%%%%%%给出在T=1时的最大误差 max( abs ( u (n+1,:) -ana (n+1,:) ) ) %图示化 [tt,xx]= meshgrid(t,x); xx=xx'; tt=tt'; mesh(tt,xx,ana);title('analytical solution'); xlabel(sprintf('steps in time domain\nfrom 0 to 1'));ylabel(sprintf('steps in x-domain\nfrom 0 to 1')); zlabel('U(x,y)'); pause mesh(tt,xx,u); title('numerical solution'); xlabel(sprintf('steps in time domain\nfrom 0 to 1'));ylabel(sprintf('steps in x-domain\nfrom 0 to 1')); zlabel('u(xi,tk)'); pause mesh(tt,xx,u-ana); title('diffeence beween analical and numeical soluion'); xlabel(sprintf('steps intime- domain\nfrom 0 to 1')); ylabel(sprintf('steps in x domain\nfrom 0 to 1')); zlabel('u-ana'); pause plot(x,u(n+1,:),'-or',x,ana(n+1,:)) title('显格式数值解与精确解的比较') legend('uh(x,1)','u(x,1)') function pdegui(action) %PDEGUI Demonstrate soluton of model partial differential equations. % PDEGUI demonstrates the finite difference solution of model problems % involving Laplace's operator: % del^2_U =partial^2_U/partial_x^2 + partial^2_U/partial_y^2. % The PDEs are: % Poisson equation, del^2_U = 1. % Eigenvalue equation, del^2_U + lambda U = 0. % Heat equation, partial_U/dt = del^2_U. % Wave equation, partial^2_U/partial_t^2 =del^2_U. % Various regions can be selected: % Square, -1 = x = 1, -1 = y = 1. % L-shaped MathWorks logo made from 3/4 of the square. % Curved 'L', with a quarter circle in the 4-th square. % Disc of radius one. % Annulus with outer radius = 1, inner radius = 1/3; % Heart-shaped cardioid. % Exterior of a "Butterfly". % Cross section of a double ridge microwave waveguide. % For the eigenvalue problem, you can set the eigenvalue index. % For the heat equation, you can set sigma = (delta_t)/(delta_x)^2. % For the wave equation, you can set sigma = ((delta_t)/(delta_x))^2. % If sigma is too large, the time-stepping methods are unstable. if nargin == 0 Initialize_GUI action ='restart'; end set(findobj('tag','stop'),'userdata',0) drawnow if isequal(action,'restart') % Retrieve parameters from GUI r = findobj('tag','region'); s = get(r,'string'); region = s{get(r,'value')}; eqn = get(findobj('tag','eqn'),'value'); n =eval(get(findobj('tag','edit1'),'string')); text2 = findobj('tag','text2'); edit2 =findobj('tag','edit2'); if eqn == 1 set(text2,'vis','off'); set(edit2,'vis','off');set(findobj('string','+'),'vis','off'); set(findobj('string','-'),'vis','off'); else set(text2,'vis','on'); set(edit2,'vis','on'); set(findobj('string','+'),'vis','on'); set(findobj('string','-'),'vis','on'); end sigma = []; indx = []; switch eqn case 1 case 2 set(text2,'string','index =')set(edit2,'string','1') indx = 1; case 3 set(text2,'string','sigma =') set(edit2,'string','0.250') sigma = .250; case 4 set(text2,'string','sigma =') set(edit2,'string','0.500') sigma = .500; end % Generate the grid. G = gengrid(region(1),n); % Generate the discrete Laplacian.A = Laplacian(G); m = size(A,1); % Initialize the solution [x,y] = meshgrid(-1:2/n:1); if eqn == 1; u = []; v = []; lambda = []; elseif eqn == 2 v = []; h = 2/n; [u,lambda] = eigswatch(-(1/h^2)*A,10); else k = (G > 0); u = zeros(m,1); v = zeros(m,1); r =sqrt((x(k)-.5).^2 + (y(k)+.5).^2); if eqn == 3 u(:) = -(15/n)*exp(-5*r); else v(:) = -(15/n)*exp(-5*r); end lambda = []; end % Save everything in figure's userdata. data.eqn = eqn; data.G = G; data.A = A; data.x = x; data.y = y; data.n = n; data.m = m; data.u = u; data.v = v; data.sigma = sigma; data.indx = indx; mbda = lambda;set(gcf,'userdata',data); else % New value of sigma or index. data = get(gcf,'userdata'); edit2 = findobj('tag','edit2'); if data.eqn == 2 data.indx = eval(get(edit2,'string')); ifisequal(action,'+') data.indx = data.indx + 1; elseif isequal(action,'-') data.indx =max(1,data.indx - 1); end set(edit2,'string',sprintf('%6.0f',data.indx)) else data.sigma = eval(get(edit2,'string')); if isequal(action,'+') data.sigma = data.sigma + .001; elseif isequal(action,'-') data.sigma = data.sigma - .001; endset(edit2,'string',sprintf('%6.3f',data.sigma)) end set(gcf,'userdata',data); end data =get(gcf,'userdata'); A = data.A; u = data.u; v = data.v; sigma = data.sigma; eqn = data.eqn; m = data.m; n = data.n; indx = data.indx; x = data.x; y = data.y; G = data.G; set(findobj('tag','stop'),'userdata',0); while get(findobj('tag','stop'),'userdata') == 0 switch eqn case 1 % Solve sparse linear system Au = 1 u = A\(ones(m,1)); u = u/max(max(abs(u))); set(findobj('tag','stop'),'userdata',1) case 2 % Some eigenvalues already computed; maybe need more. while indx >length(mbda) k = length(mbda)+10; h = 2/n; [u,lambda] = eigswatch(-(1/h^2)*A,k); mbda = lambda; data.u = u; end lambda = mbda(indx); u = data.u(:,indx); u = u/max(max(abs(u))); set(findobj('tag','stop'),'userdata',1) case 3 % One step for heat equation u = u + sigma*A*u; data.u = u; case 4 % One step for wave equation w = u; u = 2*u - v + sigma*A*u; v = w; data.u = u; data.v = v; end % Map the solution onto the grid and display it. z = G; k = find(G>0); z(k) = u(z(k));contour(x,y,z,-1:.05:1) axis square set(gca,'color','k') if eqn == 2title(sprintf('lambda(%2d) = %9.4f',indx,lambda)) end drawnow set(gcf,'userdata',data) end % ------------------------------------------------------------------------ functionInitialize_GUI clf reset set(gcf,'pos',[150 260 710420],'doublebuffer','on','menubar','none', ... 'numbertitle','off','name','PDEgui','userdata',1); set(gca,'pos',[.11 .11 .6 .815]) uicontrol( ... 'tag','region', ...'style','list', ... 'units','normal', ... 'position',[.75 .55 .20 .35], ... 'string',{'Square','L (Logo)','Curved L','Disc','Annulus',... 'Heart','Butterfly','Waveguide'}, ... 'fontsize',12, ... 'value',2, ... 'callback','pdegui(''restart'')') uicontrol( ... 'tag','eqn', ... 'style','list', ...'units','normal', ... 'position',[.75 .35 .20 .18], ...'string',{'Poisson','Eigenvalue','Heat','Wave'}, ... 'fontsize',12, ... 'value',1, ...'callback','pdegui(''restart'')') uicontrol( ... 'tag','text1', ... 'style','text', ... 'units','normal', ... 'position',[.75 .27 .10 .06], ... 'string','grid size =', ... 'fontsize',12) uicontrol( ...'tag','edit1', ... 'style','edit', ... 'units','normal', ... 'position',[.85 .27 .10 .06], ...'string','32', ... 'backgroundcolor',[1 1 1], ... 'fontsize',12, ... 'callback','pdegui(''restart'')') uicontrol( ... 'tag','text2', ... 'style','text', ... 'units','normal', ... 'position',[.75 .19 .10 .06], ... 'string','', ... 'fontsize',12) uicontrol( ... 'tag','edit2', ... 'style','edit', ... 'units','normal', ... ' position',[.85 .19 .08 .06], ... 'string','0.5', ... 'backgroundcolor',[1 1 1], ... 'fontsize',12, ... 'callback','pdegui(''sigint'')') uicontrol( ... 'style','pushbutton', ... 'units','normal', ...'position',[.93 .19 .02 .03], ... 'string','-', ... 'fontsize',12, ... 'callback','pdegui(''-'')') uicontrol( ... 'style','pushbutton', ... 'units','normal', ... 'position',[.93 .22 .02 .03], ...'string','+', ... 'fontsize',12, ... 'callback','pdegui(''+'')') uicontrol( ... 'tag','stop', ...'style','pushbutton', ... 'units','normal', ... 'position',[.80 .11 .12 .06], ... 'string','stop', ... 'fontsize',12, ... 'userdata',0, ... 'callback','set(gcbo,''userdata'',1)'); uicontrol( ...'style','pushbutton', ... 'units','normal', ... 'position',[.80 .03 .12 .06], ... 'string','close', ... 'fontsize',12, ... 'callback','close(gcf)'); % ------------------------------------------------------------------------ function [u,lambda] = eigswatch(A,k) % [u,lambda] = eigswatch(A,k) % Modify pointer and text while computing k smallest % eigenvalues and eigenvectors of sparse matrix A. ps = get(gcf,'pointer'); set(gcf,'pointer','watch') text2 =findobj('tag','text2'); edit2 = findobj('tag','edit2'); str = get(text2,'string');set(text2,'string','computing') set(edit2,'vis','off') drawnow opts.disp = 0; [u,lambda] = eigs(A,k,0,opts); [lambda,p] = sort(diag(lambda)); u = u(:,p); % Make first eigenvector positive to % distinguish it from Poisson solution. u(:,1) = abs(u(:,1));set(text2,'string',str) set(edit2,'vis','on') set(gcf,'pointer',ps) % ------------------------------------------------------------------------ function G = gengrid(R,n) %GENGRID Number the grid points in a two dimensional region. % G = GENGRID('R',n) numbers the points on an n-by-n grid in % the subregion of -1=x=1 and -1=y=1 determined by 'R'. %SPY(GENGRID('R',n)) plots the points. % LAPLACIAN(GENGRID('R',n)) generates the 5-point discrete Laplacian. % The regions currently available are: % 'S' - the entire square. % 'L' - the L-shaped domain made from 3/4 of the entire square. % 'C' - like the 'L', but with a quarter circle in the 4-th square. % 'D' - the unit disc. % 'A' - an annulus. % 'H' - a heart-shaped cardioid. % 'B' - the exterior of a "Butterfly". % 'W' - cross section of a double ridge microwave waveguide. [x,y] = meshgrid(-1:2/n:1); switch(R) case 'S' G = (x > -1) & (x 1) & (y > -1) & (y 1); case 'L' G = (x > -1) & (x 1) & (y > -1) & (y 1) & ( (x > 0) | (y > 0)); case 'C' G = (x > -1) & (x 1) & (y > -1) & (y 1) & ((x+1).^2+(y+1).^2 > 1); case 'D' G = x.^2 + y.^2 1; case 'A' G = ( x.^2 + y.^2 1) & ( x.^2 + y.^2 > 1/3); case 'H' rho = .75; sigma= .75; G = (x.^2+y.^2).*(x.^2+y.^2-sigma*y) rho*x.^2; case 'B' t = atan2(y,x); r =sqrt(x.^2 + y.^2); G = (r >= sin(2*t) + .2*sin(8*t)) & ... (x > -1) & (x 1) & (y > -1) & (y 1); case 'W' G = (x > -1) & (x 1) & (y > -1) & (y 1) & (abs(x)>1/2 | abs(y)<1/2); otherwise error('Invalid region type.'); end k = find(G); G = zeros(size(G)); G(k) = (1:length(k))'; % ------------------------------------------------------------------------ function D = Laplacian(G) %LAPLACIAN Construct five-point finite difference Laplacian. % Laplacian(G) is the sparse form of the two-dimensional, % 5-point discrete Laplacian。