final2.m

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

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