Solving system of PDES using Method of Lines

I am trying to solve 2 system of pdes using method of lines, but unfortunately im not sure why i have this error "Index exceeds the number of array elements. Index must not exceed 2.

clear all; clc; close all;
L=500;                                                                  %height of the boundary layer
z=linspace(20,L,100);                                                   %domain discretization in 100 steps
t0=0;
tf=20;
nt=100;
t=linspace(t0,tf,nt);                                                   %time span
n=numel(z);                                                             %returns number of elements in z
time=t;

%Initial Conditions
Tsta = 0.15;                                                            %°C
wtheta_0 = 0.15;                                                        %°C ms^-1
h0 = 200;                                                               %initial height of BL(metres)
h = 2*h0;                                                               %height of BL
T0=(Tsta/h).*(-0.55).*((z./h).^(-4/3)).*(1-(z./h));
wtheta_0=wtheta_0*(1-(z./h));
U0=[T0; wtheta_0];                                                       %initial values as vector

%Setting up the solver
[t,U]=ode15s(@(t,U) F(t,U,z,n),[t0,tf],U0);

Unable to perform assignment because the left and right sides have a different number of elements.

Error in solution>F (line 124)
dWthetadt(i) = A5.*(d2Wthetadz2(i))-(2.*C6).*(Wtheta(i)./tau_p)-(w_sq.*T(i)).*((1-C7).*(g*alpha.*w_sq));

Error in solution (line 21)
[t,U]=ode15s(@(t,U) F(t,U,z,n),[t0,tf],U0);

Error in odearguments (line 92)
f0 = ode(t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode15s (line 153)
odearguments(odeIsFuncHandle, odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

%t is time
%U is T and Wtheta
%z is the vector of boundary layer height

%Saving results
T=U(1:n,:);
wtheta=U(n+1:2*n,:);

%plotting


%Defining Function for method of lines
function dUdt = F(t,U,z,n)

%Initializing Derivatives: Preallocation
dTdt = zeros(n,1);
dWthetadt = zeros(n,1);
dUdt=zeros(n,1);
dTdz = zeros(n,1);
dWthetadz = zeros(n,1);
d2Wthetadz2 = zeros(n,1);

%Initializing Initial Values
T(1:n) = U(1:n);
Wtheta(1:n) = U(n+1:2*n);

%Boundary Condition
h0 = 200;
h = 2*h0;
Wtheta_0 = 0.15;
Tsta = 0.15;
zb = 0.1*h0;                                                            %Lower Boundary height
S0 = 13.4;                                                              %Stratification Parameter
gamma_i = (S0*Tsta)/h0;                                                 %Inversion lapse rate
T(:,1) = (Tsta/h)*(-0.55).*((zb./h).^(-4/3))*(1-(zb./h));               %Lower Boundary condition specified at zl=0.1*h 
T(:,n+1) = gamma_i;                                                     %Upper Boundary condition

Wtheta(:,n+1) = Wtheta_0.*(1-(zb./h));                                  %Lower Boundary condition 
Wtheta(:,2*n+1) = 0;                                                    %Upper Boundary condition

for i=2:n
    dz = 100./n;
    dTdz(i) = 1./(2.*dz).*(T(i+1)-T(i-1));                              %Centered difference
    dWthetadz(i) = 1./(2.*dz).*(Wtheta(i+1)-Wtheta(i-1));               %Centered difference
    d2Wthetadz2(i) = 1./(dz.^2).*(Wtheta(i+1)-2.*Wtheta(i)+Wtheta(i-1));%Centered difference
end
for i=2:n
    %parameters

h0=200;
h=2*h0;
B1=19.3;
l0=0.15*h;
k=0.4; %von karma constant
lamdha=0.225;
g=9.81;   %9.81m/s^2
alpha=214e-6;   %0.000214°c
zeta=-5*z/h;
tau=B1;     %Timescale
tau_p = tau;


%Derived Constants
a8=1.4286e-1;
a9=2.3810e-2;

%Constants  
C2=1;   
C4=1.75;
C5=0.3;
C6=3.25;
C7=0.5;

theta_sq=-(2.*Wtheta(i).*T(i))*tau/4*C2;
E=2*(C2/tau).*theta_sq;                                                 %dissipation rate
w_sq=-(tau_p/C4)*(2-(4/3)*C5)*g*alpha.*Wtheta(i)-4*E/3;

%q=sqrt(3*w_sq);


%surface length
if zeta >= 1            
    lsurf=(k*z)/3.7;
elseif 0 <= zeta & zeta < 1 % SH inequality reversed
    lsurf=k*z*(1+2.7*zeta)^(-1);
elseif zeta < 0
    lsurf=k.*z.*(1-100*zeta).^0.2;
end

N2 = g*alpha.*T(i);                                                    
if N2 <= 0
    lB = 1e14; Cw = 0;
else
    lB = 0.53*sqrt(3*w_sq)/sqrt(N2); Cw = 0.04;
end


l=1./((1./l0)+(1./lsurf)+(1./lB));  % SH : Brackets missing


A5=(a8.*w_sq) + (a9.*lamdha.*Wtheta(i))*tau./l;
    dTdt(i) = -d2Wthetadz2(i);                                          %Equations T and Wtheta
    dWthetadt(i) = A5.*(d2Wthetadz2(i))-(2.*C6).*(Wtheta(i)./tau_p)-(w_sq.*T(i)).*((1-C7).*(g*alpha.*w_sq));
end
dUdt = [dTdt; dWthetadt];
end

"

[Torsten]

Your vector of initial values U0 must be of size (200x1), not (2x1).

[Nuel]

i've now made changes to the IC but then but unfortunately resulted in another error" Unable to perform assignment because the left and right sides have a different number of elements."

[Torsten]

z is an array.

Thus lsurf is an array.

Thus l is an array.

Thus A5 is an array.

Thus A5.*(d2Wthetadz2(i))-(2.*C6).*(Wtheta(i)./tau_p)-(w_sq.*T(i)).*((1-C7).*(g*alpha.*w_sq)) is an array.

But dWthetadt(i) is a scalar.

So you assign a vector to a scalar.

This gives an error.

[Nuel]

i'm not sure exactly where ive gone wrong, but dwthetadt should be a vector

[Torsten]

i'm not sure exactly where ive gone wrong, but dwthetadt should be a vector

Yes, but not its i-th element. As I deduced above, the problem starts when you define lsurf using the complete z-vector. Maybe you only have to use z(i) instead of z in the statements

zeta=-5*z/h;

and

if zeta >= 1            
    lsurf=(k*z)/3.7;
elseif 0 <= zeta & zeta < 1 % SH inequality reversed
    lsurf=k*z*(1+2.7*zeta)^(-1);
elseif zeta < 0
    lsurf=k.*z.*(1-100*zeta).^0.2;
end

Answer 1

[Torsten]

Try this code. At least, it gives a result.

clear all; clc; close all;
L=500;                                                                  %height of the boundary layer
z=linspace(20,L,100);                                                   %domain discretization in 100 steps
t0=0;
tf=20;
nt=100;
t=linspace(t0,tf,nt);                                                   %time span
n=numel(z);                                                             %returns number of elements in z
time=t;

%Initial Conditions
Tsta = 0.15;                                                            %°C
wtheta_0 = 0.15;                                                        %°C ms^-1
h0 = 200;                                                               %initial height of BL(metres)
h = 2*h0;                                                               %height of BL
T0=(Tsta/h).*(-0.55).*((z./h).^(-4/3)).*(1-(z./h));
wtheta_0=wtheta_0*(1-(z./h));


U0=[T0, wtheta_0].';                                                       %initial values as vector

%Setting up the solver
[t,U]=ode15s(@(t,U) F(t,U,z,n),[t0,tf],U0);
U = U.';

%t is time
%U is T and Wtheta
%z is the vector of boundary layer height

%Saving results
T=U(1:n,:);
wtheta=U(n+1:2*n,:);
figure(1)
plot(z.',T)

figure(2)
plot(z.',wtheta)

%plotting


%Defining Function for method of lines
function dUdt = F(t,U,z,n)

%Initializing Derivatives: Preallocation
dTdt = zeros(n,1);
dWthetadt = zeros(n,1);
dUdt=zeros(n,1);
dTdz = zeros(n,1);
dWthetadz = zeros(n,1);
d2Wthetadz2 = zeros(n,1);

%Initializing Initial Values
T(1:n) = U(1:n);
Wtheta(1:n) = U(n+1:2*n);

%Boundary Condition
h0 = 200;
h = 2*h0;
Wtheta_0 = 0.15;
Tsta = 0.15;
zb = 0.1*h0;                                                            %Lower Boundary height
S0 = 13.4;                                                              %Stratification Parameter
gamma_i = (S0*Tsta)/h0;                                                 %Inversion lapse rate
T(:,1) = (Tsta/h)*(-0.55).*((zb./h).^(-4/3))*(1-(zb./h));               %Lower Boundary condition specified at zl=0.1*h 
T(:,n+1) = gamma_i;                                                     %Upper Boundary condition

Wtheta(:,n+1) = Wtheta_0.*(1-(zb./h));                                  %Lower Boundary condition 
Wtheta(:,2*n+1) = 0;                                                    %Upper Boundary condition

for i=2:n
    dz = 100./n;
    dTdz(i) = 1./(2.*dz).*(T(i+1)-T(i-1));                              %Centered difference
    dWthetadz(i) = 1./(2.*dz).*(Wtheta(i+1)-Wtheta(i-1));               %Centered difference
    d2Wthetadz2(i) = 1./(dz.^2).*(Wtheta(i+1)-2.*Wtheta(i)+Wtheta(i-1));%Centered difference
end
for i=2:n
    %parameters

h0=200;
h=2*h0;
B1=19.3;
l0=0.15*h;
k=0.4; %von karma constant
lamdha=0.225;
g=9.81;   %9.81m/s^2
alpha=214e-6;   %0.000214°c
zeta=-5*z(i)/h;
tau=B1;     %Timescale
tau_p = tau;


%Derived Constants
a8=1.4286e-1;
a9=2.3810e-2;

%Constants  
C2=1;   
C4=1.75;
C5=0.3;
C6=3.25;
C7=0.5;

theta_sq=-(2.*Wtheta(i).*T(i))*tau/4*C2;
E=2*(C2/tau).*theta_sq;                                                 %dissipation rate
w_sq=-(tau_p/C4)*(2-(4/3)*C5)*g*alpha.*Wtheta(i)-4*E/3;

%q=sqrt(3*w_sq);


%surface length
if zeta >= 1            
    lsurf=(k*z(i))/3.7;
elseif 0 <= zeta & zeta < 1 % SH inequality reversed
    lsurf=k*z(i)*(1+2.7*zeta)^(-1);
elseif zeta < 0
    lsurf=k.*z(i).*(1-100*zeta).^0.2;
end

N2 = g*alpha.*T(i);                                                    
if N2 <= 0
    lB = 1e14; Cw = 0;
else
    lB = 0.53*sqrt(3*w_sq)/sqrt(N2); Cw = 0.04;
end


l=1./((1./l0)+(1./lsurf)+(1./lB));  % SH : Brackets missing


A5=(a8.*w_sq) + (a9.*lamdha.*Wtheta(i))*tau./l;
    dTdt(i) = -d2Wthetadz2(i);                                          %Equations T and Wtheta
    dWthetadt(i) = A5.*(d2Wthetadz2(i))-(2.*C6).*(Wtheta(i)./tau_p)-(w_sq.*T(i)).*((1-C7).*(g*alpha.*w_sq));
end
dUdt = [dTdt; dWthetadt];
end

[Nuel]

Hi Torsten, Thank you once again.

I've been trying to specify the Initial condition over some range of values of z, but unfortunately, i again ran into an error message "Unable to perform assignment because the left and right sides have a different number of elements.".

for i=1:n

if (z(i)<h0)

T0(i)=(Tsta/h).*(-0.55).*((z(i)./h).^(-4/3)).*(1-z(i)./h)).';

Wtheta0(i)=Wtheta_0.*(1-z(i)./h)).';

else

T0(i)=0.';

Wtheta0(i)=0.';

end

U0(i)=[T0,Wtheta0].';

end

[Torsten]

T0 = zeros(1,n);
Wtheta0 = zeros(1,n);
for i=1:n
    if z(i)<h0
        T0(i)=(Tsta/h).*(-0.55).*((z(i)./h).^(-4/3)).*(1-z(i)./h));
        Wtheta0(i)=Wtheta_0.*(1-z(i)./h));
    end
end
U0 = [T0,Wtheta0].';

[Nuel]

Thank You.