i want the area under a curve which varies in time.

The code is:

function ML_Diff

clc
clear all
%This is the fucking file of the exam exercise
%Explicit Euler Method is used for solving the following
%initial-boundary-value problem:
%Ut+v1*Ux=alpha1Uxx     Ut+v2*Ux=alpha2Uxx
%U(z,0)=0                                                                      %Cond iniziale
%U(0,t)=g1(t)       U(L,t)=g2(t)                                               %Cond. Al bordo
%Problem parameters
nt = 100;                                                                      %numero punti griglia asse x
nz = 150;                                                                      %numero punti griglia asse y
tc = 1;                                                                        %è il tempo caratteristico
T = 1000;                                                                       %è il tempo adimensionalizzato t/tc 
Z = 50;                                                                        %profondità di tutto il multistrato
dZ = Z/nz;                                                                     %è l'intervallo spaziale tra 2 punti zi
ni=1*10^-6;
%---------------------layer 1----------------------
z1 = Z/2;                                                                      %z1 è la z barra del layer 1 ovvero z1=(L/2)/L
alpha1= 0.05;%0.05
v1 = 0.050;%0.050;
Pe1 = (v1*z1)/alpha1;
Dm1=(z1^2)/Pe1;
Re1=(v1*z1)/ni;
%---------------------Interface--------------------
zh=Z/2;                                                                        %profondità dell'interfaccia.
%---------------------layer 2----------------------
z2 = Z/2;                                                                      %z1 è la z barra del layer 1 ovvero z1=(L/2)/L
alpha2=0.1;
v2 = 0.050;%0.005;
Pe2 = (v2*z2)/alpha2;
Dm2=(z2^2)/Pe2;
Re2=(v2*z2)/ni;
%Check of problem parameters
if any([alpha1 alpha2 nz-2 nt-2 z1 z2]<=0)
    error('Check problem parameters')
end
%Stability 
st = 1;
while (2*max(alpha1,alpha2)*(T/nt/dZ^2)+(T/nt/dZ)*abs(max(v1,v2)))>st
    nt=nt+1;
end
%Post-stabilità andiamo allora ad esplicitare i seguenti termini:
dT=T/nt;                                                                        % Non abbiamo inserito anche il delta_Z perchè questo non
%cambia, cioè è lo stesso definito ad inizio codice.
r=(dT/dZ^2); s=(dT/dZ); 
%Inizialization
h=(zh/dZ)+1;  h=double(uint16(abs(h)));                                         %index of interface
z=linspace(0,Z,nz+1); t=linspace(0,T,nt+1);
zh=z(h);
u=getPisciazz(z,nz,Z); U=u;                                                     %Fissiamo la nostra funzione concentrazione(u) pari prorpio alla
p1=alpha1*r; 
d1=v1*s;
p2=alpha2*r; 
d2=v2*s;
%Subfunction [Per le condizioni iniziali]
function Pisciazz = getPisciazz(z,nz,Z)
Pisciazz(1:nz+1,1)=0;
%Pisciazz(1:h,1)=0                                                               %Cond iniziale asse z, 1°Strato
%Pisciazz(h+1:nz+1,1)=0                                                          %Cond iniziale asse z, 2°Strato
end
%Subfunction [Per bordo z=0] 
function VarTempo = getVarTempo(t)
%VarTempo=2
Dt=100; C0=50;
if t<=Dt
    VarTempo= C0;
else
    VarTempo= 0;
end
end
%Subfunction [Per bordo z=L]
function NoemannNothing = getNoemannNothing(t)
NoemannNothing= u(end-1);
%u(end-u)
%noemann nulla du/dz=0 u(end)=u(end-u)
%end
%Finite Difference Method [Forward for the time and Central for the space]
for j=2:nt+1
    u(2:h-1)=u(2:h-1)-d1*(u(3:h)-u(2:h-1))+p1*(u(3:h)-2*u(2:h-1)+u(1:h-2));
    u(h+1:end-1)=u(h+1:end-1)-d2*(u(h+2:end)-u(h+1:end-1))+p2*(u(h+2:end)-2*u(h+1:end-1)+u(h:end-2));
    u(h)=(Dm1*(-4*u(h-1)+u(h-2))-(Dm2*(4*u(h+1)-u(h+2))))/(v1*2*dZ-v2*2*dZ-3*Dm1-3*Dm2);
    %u(h)=(((4*u(h-1)-u(h-2)))*-Dm1-(4*u(h+1)-u(h+2))*Dm2)/((v1-v2)*2*dZ-3*(Dm1+Dm2));
    u(1)=getVarTempo(t(j));
    u(end)=getNoemannNothing(t(j));
    plot(u,z,'r*:',U,z,'k');
    set(gca,'YDir','reverse');
    xlabel('C(z,t)'); ylabel('z');
    title(['t=',num2str(t(j)),'seconds']);
    axis([0 35 0 50]);
    
    pause(0.01);
end
legend(['t=',num2str(T)],'t=0');
end

Please help me faster! I have to send the project to my teacher!

[Torsten]

You want the area under a curve ? Which curve ? Which area ?

[Alfonso]

Of the curve (u,Z) as time varies, for example:

[Torsten]

u_integrated(j) = trapz(z,u)  

[Alfonso]

Torsten how can i also plot in the title this value??

If i can ask you another question, in the code there i a problem thanks to this the area is not costant (but it should be for the principle of conservation of mass). Can you help me to find the error?

[Torsten]

Maybe I could help if you arrange the code above to work (at least technically) in MATLAB online.

[Alfonso]

how can i do it?

[Torsten]

Click on "edit" , make appropriate changes (because as is, your code doesn't work ; at least it does not produce any plots) and then click on the green arrow (Run).

The code must be arranged as

script part

function1

function2

...

function n

All functions have to be closed by an "end" statement.

Nested functions are not allowed.

[Alfonso]

clc
clear all

%This is the fucking file of the exam exercise
%Explicit Euler Method is used for solving the following
%initial-boundary-value problem:
%Ut+v1*Ux=alpha1Uxx     Ut+v2*Ux=alpha2Uxx
%U(z,0)=0                                                                      %Cond iniziale
%U(0,t)=g1(t)       U(L,t)=g2(t)                                               %Cond. Al bordo

%Problem parameters
nt = 100;                                                                      %numero punti griglia asse x
nz = 150;                                                                      %numero punti griglia asse y
tc = 1;                                                                        %è il tempo caratteristico
T = 1000;                                                                       %è il tempo adimensionalizzato t/tc
Z = 50;                                                                        %profondità di tutto il multistrato
dZ = Z/nz;                                                                     %è l'intervallo spaziale tra 2 punti zi
ni=1*10^-6
%---------------------layer 1----------------------
z1 = Z/2                                                                       %z1 è la z barra del layer 1 ovvero z1=(L/2)/L
alpha1= 0.05;%0.05
v1 = 0.050%0.050;
Pe1 = (v1*z1)/alpha1;
Dm1=(z1^2)/Pe1;
Re1=(v1*z1)/ni
%---------------------Interface--------------------
zh=Z/2;                                                                        %profondità dell'interfaccia.
%---------------------layer 2----------------------
z2 = Z/2                                                                       %z1 è la z barra del layer 1 ovvero z1=(L/2)/L
alpha2=0.1;
v2 = 0.050%0.005;
Pe2 = (v2*z2)/alpha2;
Dm2=(z2^2)/Pe2
Re2=(v2*z2)/ni
%Check of problem parameters
if any([alpha1 alpha2 nz-2 nt-2 z1 z2]<=0)
    error('Check problem parameters')
end

%Stability
st = 1;
while (2*max(alpha1,alpha2)*(T/nt/dZ^2)+(T/nt/dZ)*abs(max(v1,v2)))>st
    nt=nt+1;
end

%Post-stabilità andiamo allora ad esplicitare i seguenti termini:
dT=T/nt;                                                                        % Non abbiamo inserito anche il delta_Z perchè questo non
                                                                                %cambia, cioè è lo stesso definito ad inizio codice.
r=(dT/dZ^2); s=(dT/dZ);

%Inizialization
h=(zh/dZ)+1;  h=double(uint16(abs(h)));                                         %index of interface
z=linspace(0,Z,nz+1); t=linspace(0,T,nt+1);
zh=z(h);
u=getPisciazz(z,nz,Z); U=u;                                                     %Fissiamo la nostra funzione concentrazione(u) pari prorpio alla
A(1) = trapz(z,u)
p1=alpha1*r;
d1=v1*s;
p2=alpha2*r;
d2=v2*s;



%Finite Difference Method [Forward for the time and Central for the space]
for j=2:nt+1
    u(2:h-1)=u(2:h-1)-d1*(u(3:h)-u(2:h-1))+p1*(u(3:h)-2*u(2:h-1)+u(1:h-2));
    u(h+1:end-1)=u(h+1:end-1)-d2*(u(h+2:end)-u(h+1:end-1))+p2*(u(h+2:end)-2*u(h+1:end-1)+u(h:end-2));
    u(h)=(Dm1*(-4*u(h-1)+u(h-2))-(Dm2*(4*u(h+1)-u(h+2))))/(v1*2*dZ-v2*2*dZ-3*Dm1-3*Dm2);
    %u(h)=(((4*u(h-1)-u(h-2)))*-Dm1-(4*u(h+1)-u(h+2))*Dm2)/((v1-v2)*2*dZ-3*(Dm1+Dm2));
    u(1)=getVarTempo(t(j));
    u(end)=getNoemannNothing(t(j),u);
    plot(u,z,'r*:',U,z,'k');
    set(gca,'YDir','reverse');
    xlabel('C(z,t)'); ylabel('z');
    title(['t=',num2str(t(j)),'seconds']);
    axis([0 35 0 50]);
    pause(0.01);
    A(j) = trapz(z,u);

end

legend(['t=',num2str(T)],'t=0');
figure(2)
plot(0:dT:T,A)

%Subfunction [Per le condizioni iniziali]
    function Pisciazz = getPisciazz(z,nz,Z)
        Pisciazz(1:nz+1,1)=0;
        %Pisciazz(1:h,1)=0                                                               %Cond iniziale asse z, 1°Strato
        %Pisciazz(h+1:nz+1,1)=0                                                          %Cond iniziale asse z, 2°Strato
    end
%Subfunction [Per bordo z=0]
function VarTempo = getVarTempo(t)
%VarTempo=2
Dt=100; C0=50;
if t<=Dt
    VarTempo= C0;
else
    VarTempo= 0;
end
end
%Subfunction [Per bordo z=L]
function NoemannNothing = getNoemannNothing(t,u)
NoemannNothing= u(end-1);
%u(end-u)
%noemann nulla du/dz=0 u(end)=u(end-u)
end

[Torsten]

I modified the code to include the area computation (see above).

Since I can't decode what you have programmed as PDE, initial and boundary conditions, I cannot help you regarding possible errors.

[Alfonso]

The code should show the course of a contaminant in a multilayer soil. It is assumed that the contaminant is introduced for a time of Dt:=100 seconds, I have then used an impulsive Dirichlet as a condition at the upper edge (on the ground plane). Then at the interface between the 2 layers there is an equilibrium condition of contaminant flow and finally at the lower edge (Z=50m) there is a null Normanna condition. The equation used is the advection-diffusion equation for a usual soil dispersion and the approximations used are: *Forward for space (derivative 1⁰) *Forward for time *Central for the 2⁰ derivative of space

Mass should be conserved but unfortunately I have not understood why mass is not conserved

[Torsten]

And what do you expect as result for the array A = trapz(z,u) from mass conservation ?

At least it cannot be constant because you have an inflow of contaminant for 100 seconds, as you state above. Thus A must be increasing during this time period.

[Alfonso]

The fact Is that the area under the curve Is Just the Mass of contaminant, i've ran the code and discovered that the area Is not the same during che process above, but Is should be.

I think It Is a provlem for the diffusivity approximation terme, so maybe its wrong the centrale approximation, dont you think Is Is?

[Torsten]

Didn't you read my answer ? The mass of the contaminant changes with time (increases) since you inject more and more of it over a period of 100 s. So how should it be the same during the process ?