transport1d_upwind.m

来自「１Ｄ非连续性力学分析」· M 代码 · 共 52 行

%*************************
%  transport 1D upwind
%************************

clear all
close all

%*** declaration des variables

L=100; % longueur de la poutre
v=0.02; % vitesse constante

I=50; % nombre element en espace 0<i<I
dx=L/(I); % increment en espace

N=200; % temps total 0<n<N
dt=0.1*dx/v; %increment de temps



%********* condition initiale et condition limite****
phi=sparse(zeros(I,N));
for i=1:I
    x(i)=(i-1)*dx;
    phi(i,1)=1-x(i);
end
for n=2:N
    phi(1,n)=1-x(1)+v*n*dt;
end



%boucle sur temps n>2
for n=2:N
    for i=2:I
        phi(i,n)=-(v*(dt/dx)*(phi(i,n-1)-phi(i-1,n-1)))+phi(i,n-1);
        phiex(i,n)=1-x(i)+v*n*dt;
    end
end



%recherche du front stabilite et diffusion
for n=1:N
    for i=2:I
        if (phi(i,n)<0)
            k=i;
            a=(phi(k-1,n)-phi(k,n))/(x(k-1)-x(k));
            b=phi(k,n)-a*x(k);
            xk(n)=-(b/a);
            x_exacte(n)=1+v*n*dt;
            %stabilit