transport1d_leap_frog.m

１Ｄ非连续性力学分析

%**************************
%  transport 1D leap frog
%**************************

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-1); % increment en espace

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



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



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


%recherche du front stabilite et diffusion
for n=3:N
    indices=[];
    for i=2:(I-1)
        if (phi(i,n)<0)
            indices=[indices i];
        end
    end
    indice1=indices(1);
    k=indice1;
    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;
    erreur_x(n)=100*(x_exacte(n)-xk(n))/x_exacte(n);
    erreur_dif(n)=1-xk(n)+v*n*dt;
    tmp=n;  
end


%sauvegarde des resultats de stabilit