wenping.m

来自「非线性偏微分方程用不同差分格式求解两种TVD格式」· M 代码 · 共 50 行

clear;
lambda=0.5;
h=0.001;
tao=lambda*h;
t=1;
s=0.5;
x=0:h:2;
for i=1:length(x)
    if x(i)<0.6 & x(i)>0.4
        u10(i)=1;
    else
        u10(i)=0;
    end
end
u=zeros(t/tao,size(u10,2));
u(1,:)=u10;
for n=1:t/tao;
    v=zeros(size(u)); 
    for j=2:length(x)
        if u(n,j)==u(n,j-1)
            v(n,j)=u(n,j);
        elseif u(n,j-1)>0 && (u(n,j-1)+u(n,j))>0
            v(n,j)=u(n,j-1);
        elseif u(n,j)<0 && (u(n,j-1)+u(n,j))<0
            v(n,j)=u(n,j);
        else
            v(n,j)=0;
        end
    end
    for j=1:length(x)-1
        u(n+1,j)=u(n,j)-lambda*(v(n,j+1)^2-v(n,j)^2)/2;
    end
end
n=t/tao;
for i=1:length(x)
    if x(i)<0.4
        ustar(i)=0;
    elseif x(i)>0.6+s*n*tao && n*tao<=0.4 || (n*tao>=0.4 && x(i)>(sqrt(2*0.2*n*tao)+0.4))
        ustar(i)=0;
    elseif x(i)>=0.4 && x(i)<=0.4+n*tao
        ustar(i)=(x(i)-0.4)/n/tao;
    else
        ustar(i)=1;
    end
end
%plot(x,u(t/tao,:))
%hold on
plot(x,ustar,'r')
%hold off