Back

IV-MATLAB PROGRAMS

1-Program ploting V'(u):

nu=-2/27;

alpha=1;

beta=-1;

K=0.06;

u=-0.4:0.001:1;

f=(nu+alpha*u.^2+beta*u.^3)/K;

plot(u,f)

grid

2-Programs ploting the equilibriums ue(x) of our equation:

% function which return the values of u(x)
function U=Equi(nu,alpha,beta,K,dx,L,U0,D)

ktermes=L/dx;

% limits conditions u(0)=u0 and P(0)=D with the introduction of a virtual point u-1in order to conserve the stability of the scheme.

U=[U0];
U1=U0+D*dx-(dx.^2).*(nu+alpha.*U0.^2+beta.*U0.^3)./(2*K);
U=[U U1];

Uk=U1;
Ukmoins1=U0;

for k=1:ktermes-1
Ukplus1=2.*Uk-Ukmoins1 -(dx.^2).*(nu+alpha.*Uk.^2+beta.*Uk.^3)./K;
Ukmoins1=Uk;
Uk=Ukplus1;
U=[U Ukplus1];
end

% function which return the values of P
function V=DerivEqui(nu,alpha,beta,K,dx,L,U0,D)

ktermes=L/dx;

U1=U0+D*dx-(dx.^2).*(nu+alpha.*U0.^2+beta.*U0.^3)./(2*K);
V=[D];

Uk=U1;
Ukmoins1=U0;

for k=1:ktermes
Ukplus1=2.*Uk-Ukmoins1 -(dx.^2).*(nu+alpha.*Uk.^2+beta.*Uk.^3)./K;
dUk=(Ukplus1-Uk)./dx;
Ukmoins1=Uk;
Uk=Ukplus1;
V=[V dUk];
end

% program ploting the equilibriums in the plane (u,P) for differents values of D=P(0).
clf
alpha=1;
beta=1;
dx=5*10^-2;
L=15.0;
U0=-1/3+(1/3)^(1/2);
nu=-2/27;
K=0.6;
x=0:dx:L;

for D=0:0.05:0.7
U=Equi(nu,alpha,beta,K,dx,L,U0,D);
V=DerivEqui(nu,alpha,beta,K,dx,L,U0,D);
hold on
plot(U,V)
end

3-Program ploting the curves u(x,t) of our equation:

% program calculating u(x,t) according to the differents parameters of our equation and the initials and limits conditions which correspond to a profile with a gap between u(0) and u(L).
function U=cassure(nu,alpha,beta,K,dt,dx,t,L,Ue0,UeL)

U=[Ue0];

ntermes=t/dt;

ktermes=L/dx;

V=[Ue0];

for k=1:ktermes/2

Ukn=Ue0;

U=[U Ukn];

end

for k=1:ktermes/2

Ukn=UeL;

U=[U Ukn];

end

for n=1:ntermes

V=[Ue0];Ukmoins1n=Ue0;

for k=1:ktermes-1

Ukn=U(n,k+1);

Ukplus1n=U(n,k+2);

Uknplus1=Ukn+dt.*(nu+alpha.*Ukn.^2+beta.*Ukn.^3)+K.*dt.*(Ukmoins1n-2.*Ukn+Ukplus1n)./(dx.^2);

Ukmoins1n=Ukn;

V=[V Uknplus1];

end

U=[U;V UeL];

end

% program ploting the trajectories u(x,t) according to the parameters we fixe.
clf

alpha=1;

beta=1;

dt=2*10^-2;

dx=5*10^-2;

t=2.0;

L=1.0;

Ue0=-1/3+(1/3)^(1/2);

UeL=-1/3-(1/3)^(1/2);

nu=-2/27;

y=0:dt:t;

x=0:dx:L;

[X,Y]=meshgrid(x,y);

for K=0:0.01:0.06

figure

Z=cassure(nu,alpha,beta,K,dt,dx,t,L,Ue0,UeL);

surf(X,Y,Z)

view(10,60)

end

4-:Program ploting the profile u(x) for t the last value of time:

clf

alpha=1;

beta=1;

dt=5*10^-3;

dx=5*10^-2;

K=0.24;

L=1.0;

Ue0=-1/3+(1/3)^(1/2);

UeL=-1/3-(1/3)^(1/2);

nu=-2/27;

ktermes=L/dx;

x=0:dx:L;

hold on

t=20;

Z=cassure(nu,alpha,beta,K,dt,dx,t,L,Ue0,UeL);

ntermes=t/dt;

z=Z(ntermes+1,1);

for j=2:ktermes+1

z=[z,Z(ntermes+1,j)];

end

plot(x,z)

5-Program ploting the profile u(P) for t the last value of time:

clf

alpha=1;

beta=1;

dt=5*10^-3;

dx=5*10^-2;

K=0.24;

L=1.0;

Ue0=-1/3+(1/3)^(1/2);

UeL=-1/3-(1/3)^(1/2);

nu=-2/27;

ktermes=L/dx;

x=0:dx:L;

hold on

t=20;

Z=cassure(nu,alpha,beta,K,dt,dx,t,L,Ue0,UeL);

ntermes=t/dt;

z=Z(ntermes+1,1);

dz=(Z(ntermes+1,2)-Z(ntermes+1,1))./dx;

for j=3:ktermes+1

z=[z,Z(ntermes+1,j-1)];

dz=[dz,(Z(ntermes+1,j)-Z(ntermes+1,j-1))./dx];

end

z=[z,Z(ntermes+1,ktermes+1)];

dz=[dz,(Z(ntermes+1,ktermes+1)-Z(ntermes+1,ktermes))./dx];

plot(z,dz)