scheme.m

求解激波管问题等产生的偏微分方程组的各种数值算法。

script
initial
uOld=zeros(3,lenX+1);uNew=zeros(3,lenX+1);
uOld(1,:)=rho;uOld(2,:)=m;uOld(3,:)=eng;
count=0;
aviobj=avifile('result.avi','fps',3);
for time=0:deltT:endTime
    for k=2:lenX
         cof=[0,1,0;
             (gamma-3)/2*(uOld(2,k)/uOld(1,k)).^2,-(gamma-3)*uOld(2,k)/uOld(1,k),gamma-1;
            -uOld(2,k)/uOld(1,k).^2*(gamma*uOld(3,k)-(gamma-1)*uOld(2,k).^2/uOld(1,k)), ...
            1/uOld(1,k)*(gamma*uOld(3,k)-1.5*(gamma-1)*uOld(2,k).^2/uOld(1,k)),gamma*uOld(2,k)/uOld(1,k)];
         uNew(:,k)=uOld(:,k)-R/2*cof*(uOld(:,k+1)-uOld(:,k-1))  ...
                    +R.^2/2*cof*cof*(uOld(:,k+1)-2*uOld(:,k)+uOld(:,k-1));
    end
    uNew(:,1)=uNew(:,2);uNew(:,lenX+1)=uNew(:,lenX);
    uOld=uNew;
    
    vel=uNew(2,:)./uNew(1,:);pres=(gamma-1)*(uNew(3,:)-0.5*uNew(1,:).*vel.*vel);
    plot(-2:deltX:2,uNew(1,:),'-',-2:deltX:2,vel,'-.',-2:deltX:2,pres,':');
    title(strcat('time=',num2str(time)));legend('density','velcocity','pressure');  
    axis([-2 2 -0.2 2.2]);
    count=count+1;
    if(rem(count,25)==0)
         F(count) = getframe;aviobj=addframe(aviobj,F(count));
     end     
end
   aviobj=close(aviobj);
   movie(F);  % show the movie