nonlinear_lw_lf.m

非线性偏微分方程用不同差分格式求解两种TVD格式

	clear;
	
	cfl=0.5;
    h=0.01;
	dt=cfl*h;
	tn=0.2;
	n=tn/dt;
	x1=0;
	x2=0.4;
	x3=0.6;
	x4=2.5;
	m=(x4-x1)/h;
	s=0.5;
	
    	t_ex=0;
	t_lw=0;
	t_gn=0;
	t_tvd=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
        t_ex=t_ex+dt;
        if(t_ex<=0.4)
        for i=1:1:m+1
            if(x(i)>=x2&&x(i)<=t_ex+x2)
            un1_ex(i)=(x(i)-x2)/t_ex;
            elseif(x(i)>t_ex+x2&&x(i)<=s*t_ex+x3)
            un1_ex(i)=1;
            else
            un1_ex(i)=0;
            end
        end
        else
        for i=1:1:m+1
            if(x(i)>=x2&&x(i)<=s*t_ex+x3)
            un1_ex(i)=(x(i)-x2)/t_ex;
            else
            un1_ex(i)=0;
            end
        end
	end
    un_ex=un1_ex;
    end
    
	
	%LAX-WENDROFF SCHEME
 	un_lw=un;
	un1_lw=un1;
	while t_lw<tn
		t_lw=t_lw+dt;
		tmp1=un_lw(2)^2-un_lw(1)^2;
		tmpm1=un_lw(m+1)^2-un_lw(m)^2;
		un1_lw(1)=un_lw(1)-0.25*cfl*un_lw(2)^2+0.125*cfl^2*((un_lw(1)+un_lw(2))*tmp1-un_lw(1)^3);
        	un1_lw(m+1)=un_lw(m+1)+0.25*cfl*un_lw(m)+0.125*cfl^2*(-un_lw(m+1)^3-(un_lw(m+1)+un_lw(m))*tmpm1);
        	for i=2:1:m
        	tmpi1=un_lw(i+1)^2-un_lw(i-1)^2;
        	tmpi2=un_lw(i)+un_lw(i+1);
        	tmpi3=un_lw(i+1)^2-un_lw(i)^2;
        	tmpi4=un_lw(i)+un_lw(i-1);
        	tmpi5=un_lw(i)^2-un_lw(i-1)^2;
            	un1_lw(i)=un_lw(i)-0.25*cfl*tmpi1+0.125*cfl^2*(tmpi2*tmpi3-tmpi4*tmpi5);
	    	end
    	un_lw=un1_lw;
    end
	
	%GODUNOV SCHEME
	un_gn=un;
	un1_gn=un1;
	g=zeros(m+1,1);
	while t_gn<tn
		t_gn=t_gn+dt;
		for j=1:1:m
		for i=2:1:m+1
			if(un_gn(i-1)==un_gn(i))
				g(i)=un_gn(i)^2/2;
			elseif(un_gn(i-1)>0)
				g(i)=un_gn(i-1)^2/2;
			else
				g(i)=0;
			end
		end
		g(1)=un_gn(1);
		un1_gn(j)=un_gn(j)-cfl*(g(j+1)-g(j));
		un1_gn(m+1)=un_gn(m+1)-cfl*(un_gn(m+1)-un_gn(m));
        	end
    un_gn=un1_gn;
	end
	
	%TVD SCHEME
	un_tvd=un;
	un1_tvd=un1;
	theta=zeros(m+1,1);
	fai=zeros(m+1,1);
	gn=zeros(m+1,1);
    a=0;
	while t_tvd<tn
		t_tvd=t_tvd+dt;
		for i=2:1:m+1
		    if(i~=m+1)
                if(un_tvd(i+1)~=un_tvd(i))
		        theta(i)=(un_tvd(i)-un_tvd(i-1))/(un_tvd(i+1)-un_tvd(i));
                else
                theta(i)=0;
                end
		        fai(i)=(abs(theta(i))+theta(i))/(1+abs(theta(i)));
%fai(i)=1;
                a=(un_tvd(i)+un_tvd(i+1))/2;
    		    tmp1=un_tvd(i)^2/2;
		        tmp2=un_tvd(i+1)^2/2;	
                if(un_tvd(i)==0&&un_tvd(i+1)==0)
		            gn(i)=0.5*(tmp1+tmp2)-0.5/cfl*(un_tvd(i+1)-un_tvd(i));
                else
                    gn(i)=0.5*(tmp1+tmp2)-0.5/cfl*(un_tvd(i+1)-un_tvd(i))+0.5*fai(i)*(1./cfl/a-cfl*a)*(tmp2-tmp1);
                end
            else		    
		        tmp1=un_tvd(m+1)^2/2;
		        tmp2=0;
                gn(m+1)=0.5*(tmp1+tmp2)-0.5/cfl*(-un_tvd(i));
            end
		    un1_tvd(i)=un_tvd(i)-cfl*(gn(i)-gn(i-1));
		end
		un1_tvd(1)=un_tvd(i)-cfl*gn(i);
		un_tvd=un1_tvd;
	end		    
    plot(x,un1_ex,x,un1_gn,x,un1_lw,x,un1_tvd);
    %axis([1.5,2.5,-0.5,1.5]);
	
	%save('difference scheme final.mat','n','m','un1_ex','un1_uw','un1_lf','un1_lw','un1_bw');