1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43
|
clear all; close all; a=1; T=1; h=0.0025; x=[0:h:1]; tau=0.0025; t=[0:tau:T]; r=a*tau/h; M=length(x)-1; N=length(t)-1; [T X]=meshgrid(t,x);
e=r^2*ones(M-1,1); A=spdiags([e 2*(1-e) e],[-1 0 1],M-1,M-1);
u=zeros(M+1,N+1); u(:,1)=sin(4*pi*x'); u(:,2)=sin(4*pi*x')+tau*sin(8*pi*x'); u(1,:)=0; u(end,:)=0; for n=2:N u(2:M,n+1)=A*u(2:M,n)-u(2:M,n-1); u(2,n+1)=u(2,n+1)+r^2*u(1,n); u(M,n+1)=u(M,n+1)+r^2*u(end,n); end
mesh(t,x,u), view(20,40) set(gca,'fontsize',12) xlabel('t','fontsize', 14) ylabel('x','fontsize',14) zlabel('u','fontsize',14)
ue=sin(4*pi*X).*cos(4*pi*T)+sin(8*pi*X).*sin(8*pi*T)/(8*pi);
disp('t = 1, 2, 3, 4, 5') Error=max(abs(ue(:,1/tau+1:1/tau:end)-u(:,1/tau+1:1/tau:end)))
|