final.m

关于不同的差分格式求解对流偏微分方程的比较的程序

	clear;
	
	cfl=0.5;
    h=0.01;
	dt=cfl*h;
	tn=2;
	n=tn/dt;
	x1=0;
	x2=0.4;
	x3=0.6;
	x4=2;
	m=(x4-x1)/h;
	
	x=zeros(m+1,n);
    for j=1:1:n
	for i=1:1:m+1
        x(i,j)=x1+(i-1)*h;
	end
    end
	un=zeros(m+1,1);
	un1=zeros(m+1,n);
	for i=1:1:m+1
        if (x(i)>=x2 & x(i)<=x3)
        	un(i,1)=1.0;
        else
        	un(i,1)=0;
        end
	end
	
    %exact solution
    un_ex=un;
    un1_ex=un1;
    for j=1:1:n
        for i=1:1:m+1
            tmp=x(i)-j*dt;
            if(tmp>=0.4&tmp<=0.6)
                un1_ex(i,j)=1.;
            else
                un1_ex(i,j)=0.;
            end
        end
    end
    	
	%UPWIND SCHEME
	un_uw=un;
	un1_uw=un1;
	for j=1:1:n
        un1_uw(1,j)=un_uw(1)-cfl*(un_uw(1)-0);
	for i=2:1:m+1
        un1_uw(i,j)=un_uw(i)-cfl*(un_uw(i)-un_uw(i-1));
	end
	un_uw=un1_uw(:,j);
	end
	
	%LAX-FRIEDRICHS SCHEME
	un_lf=un;
	un1_lf=un1;
	for j=1:1:n
        un1_lf(1,j)=0.5*(0+un_lf(2))-0.5*cfl*(un_lf(2)-0);
        un1_lf(m+1,j)=0.5*(un_lf(m)+0)-0.5*cfl*(0-un_lf(m));
	for i=2:1:m
        un1_lf(i,j)=0.5*(un_lf(i-1)+un_lf(i+1))-0.5*cfl*(un_lf(i+1)-un_lf(i-1));
	end
	un_lf=un1_lf(:,j);
	end
	
	%LAX-WENDROFF SCHEME
    
	un_lw=un;
	un1_lw=un1;
	for j=1:1:n
        un1_lw(1,j)=un_lw(1)-0.5*cfl*(un_lw(2)-0)+0.5*cfl^2*(un_lw(2)-2*un_lw(1)+0);
        un1_lw(m+1,j)=un_lw(m+1)-0.5*cfl*(0-un_lw(m))+0.5*cfl^2*(0-2*un_lw(m+1)+un_lw(m));
        for i=2:1:m
            un1_lw(i,j)=un_lw(i)-0.5*cfl*(un_lw(i+1)-un_lw(i-1))+0.5*cfl^2*(un_lw(i+1)-2*un_lw(i)+un_lw(i-1));
	    end
        un_lw=un1_lw(:,j);
        end
	
	%BEAM-WARMING
	un_bw=un;
	un1_bw=un1;
	for j=1:1:n
	for i=3:1:m+1
        un1_bw(i,j)=un_bw(i)-cfl*(un_bw(i)-un_bw(i-1))-0.5*cfl*(1-cfl)*(un_bw(i)-2*un_bw(i-1)+un_bw(i-2));
	end
	un_bw=un1_bw(:,j);
	end
	
    
	%subplot(2,2,1),plot(x,un,'r',x,un1_uw,'b'),axis([x1 x4 -0.5 1.5]);
	%subplot(2,2,2),plot(x,un,'r',x,un1_lf,'b'),axis([x1 x4 -0.5 1.5]);
    %subplot(2,2,3),plot(x,un,'r',x,un1_lw,'b'),axis([x1 x4 -0.5 1.5]);
    %subplot(2,2,4),plot(x,un,'r',x,un1_bw,'b'),axis([x1 x4 -0.5 1.5]);
    %plot(x,un1_ex,x,un1_uw,x,un1_lf,x,un1_lw,x,un1_bw);
   
	t=[1:1:n]';
	comet(x,un1_uw(:,t),'r');
    axis([0,2,-0.2,1.2]);
    %subplot(1,2,2),comet(x,un1_uw(:,t),'g');
	%save('difference scheme final.mat','n','m','x','un1_ex','un1_uw','un1_lf','un1_lw','un1_bw');