du-fort-frankel-difference.m

一个小程序

%
% file name: Du-Fort-Frankel-Difference.m
%


%==========================================================================
%   get space grid
clc;
clear;
h=0.1;
x=0:h:1;
%==========================================================================


%==========================================================================
%   calculate exact solution
t0=0.1*ones(size(x));
u0=exp(-pi*pi*t0).*sin(pi*x)
%==========================================================================


%==========================================================================
%   calculate numerical solution with lambda=0.1
lambda1=0.1;
tau1=h*h*lambda1;
t1=0:tau1:1;
[X1,T1]=meshgrid(x,t1);
u1=zeros(size(X1));
u1(1,:)=sin(pi*x);
u1(:,1)=zeros(size(t1));
u1(:,11)=u1(:,1);
%
v1=zeros(size(x));
for i=1:11
    if (i>1)&(i<11)
        v1(1,i)=u1(1,i)+0.5*lambda1*(u1(1,i+1)-2*u1(1,i)+u1(1,i-1));
    else
        v1(1,i)=0;
    end
end
%
p1=2*lambda1/(1+2*lambda1);
q1=(1-2*lambda1)/(1+2*lambda1);
for i=2:10
    u1(2,i)=p1*(v1(1,i+1)+v1(1,i-1))+q1*u1(1,i);
end
%
for n=3:1001
    for j=2:10
        u1(n,j)=p1*(u1(n-1,j+1)+u1(n-1,j-1))+q1*u1(n-2,j);
    end
end
%
u1(101,:)
%==========================================================================
%   calculate numerical solution with lambda=0.5
lambda2=0.5;
tau2=h*h*lambda2;
t2=0:tau2:1;
[X2,T2]=meshgrid(x,t2);
u2=zeros(size(X2));
u2(1,:)=sin(pi*x);
u2(:,1)=zeros(size(t2));
u2(:,11)=u2(:,1);
%
v2=zeros(size(x));
for i=1:11
    if (i>1)&(i<11)
        v2(1,i)=u2(1,i)+0.5*lambda2*(u2(1,i+1)-2*u2(1,i)+u2(1,i-1));
    else
        v2(1,i)=0;
    end
end
%
p2=2*lambda2/(1+2*lambda2);
q2=(1-2*lambda2)/(1+2*lambda2);
for i=2:10
    u2(2,i)=p2*(v2(1,i+1)+v2(1,i-1))+q2*u2(1,i);
end
%
for n=3:201
    for j=2:10
        u2(n,j)=p2*(u2(n-1,j+1)+u2(n-1,j-1))+q2*u2(n-2,j);
    end
end
%
u2(21,:)
%==========================================================================


%==========================================================================
%   plot for comparison
subplot(1,2,1)
plot(x,u0,x,u0,'.',x,u1(101,:),'g',x,u1(101,:),'g+')
axis([0,1,0,0.4])
grid on
title('{\bullet} : exact solution; + : numerical solution with {\lambda}=0.1')
xlabel('x')
ylabel('u(x,0.1)')
%
subplot(1,2,2)
plot(x,u0,x,u0,'.',x,u2(21,:),'g',x,u2(21,:),'g+')
axis([0,1,0,0.4])
grid on
title('{\bullet} : exact solution; + : numerical solution with {\lambda}=0.5')
xlabel('x')
ylabel('u(x,0.1)')
%==========================================================================