Solving 3 systems of pde using pdepe

I have 2 sytems of PDEs as seen in the attached file. I transformed the systems into 3 systems of pde by letting w = \partial \theta/\partial z. Such that

\theta(z,t) = u(1)

w(z,t) = u(2)

v(z,t) = u(3)

and

\partial \theta/\partial z = dudx(1)

\partial w/partial z = dudx(2)

\partial v/\partial z = dudx(3)

Below is my trial code. I have also attached the systems of pde and the trial code. I want to plot v(z,t) and \theta(z,t).

Thank you.

% BC conditions:

function [pl, ql, pr, qr] = bc(xl, ul, xr, ur, t)

pl = [0; 0; 0];

ql = [1;1;1];

pr = pl

qr = ql;

end

%initial conditions

function u0 = init(x)

D = 0.0002

thetab = 0.0001

u0 = [thetab*sin(pi/D)*z ;0; 0];

end

% Next, we will define the function

function [c, f, s] = eq1(x, t, u, dudx)

% parameters.

k1 = 6*10^(-12);

eta1 = 0.0240;

apha3 = -0.001104;

gama1 = 0.1093;

% activity parameter

xi = 0.1;

% c, f and s values

c = [gama1;apha3; 0];

f = [k1*dudx(2)+ apha3*u(3); -eta1*dudx(3); 0];

s = [0;-xi*u(2); dudx(1) - u(2)];

end

% solving systems of PDE

% our problem in cartesian coordinates

m = 0;

% Effect along the river side

% Assume the distance along the river side is x

x= linspace(0, 1,50);

t = linspace(0, 10, 40);

%Solution

sol = pdepe(m, @eq1, @init, @bc, x, t)

% Director angle

u1 = sol(:, :, 1);

% velocity flow

u2 = sol(:, :, 3);

%plot the solution

subplot(2,1,1)

surf(x, t, u1)

title('Director angle')

xlabel('x')

ylabel('t')

zlabel('u(x,t)')

view([150 25])

subplot(2,1,2)

surf(x, t, u3)

xlabel('x')

ylabel('t')

title('velocity flow')

[Torsten]

% solving systems of PDE
% our problem in cartesian coordinates
m = 0;
% Effect along the river side
% Assume the distance along the river side is x
x= linspace(0, 1,50);
t = linspace(0, 10, 40);

%Solution
sol = pdepe(m, @eq1, @init, @bc, x, t)

Error using pdepe
Spatial discretization has failed. Discretization supports only parabolic and elliptic equations, with flux term involving spatial derivative.

% Do concentration
u1 = sol(:, :, 1);

% BOD concentration
u2 = sol(:, :, 3);


%plot the solution
subplot(2,1,1)
surf(x, t, u1)    
title('Director angle')
xlabel('x')
ylabel('t')
zlabel('u(x,t)')
view([150 25])
subplot(2,1,2)
surf(x, t, u3)  
xlabel('x')
ylabel('t')
title('velocity flow')

% BC conditions:
function [pl, ql, pr, qr] = bc(xl, ul, xr, ur, t)
pl = [0; 0; 0];
ql = [1;1;1];
pr = pl;
qr = ql;
end


%initial conditions
function u0 = init(x) 
D = 0.0002;
thetab = 0.0001;
u0 = [thetab*sin(pi/D)*x ;0; 0];
end

% Next, we will define the function, init cond. and BCs.
function [c, f, s] = eq1(x, t, u, dudx)

% parameters.
k1 = 6*10^(-12);
eta1 = 0.0240;
apha3 = -0.001104;
gama1 = 0.1093;

% activity parameter
xi = 0.1; 

% c, f and s values
c = [gama1;apha3; 0];
f = [k1*u(2)+ apha3*u(3); -eta1*dudx(3); u(1)];
s = [0;-xi*u(2); u(2)];
end

[University Glasgow]

Please how can I fix the error. I'm using matlab for the first time.

Error using pdepe

Spatial discretization has failed. Discretization supports only parabolic and elliptic equations, with flux term involving spatial derivative.

[Torsten]

Your system of PDEs does not fit in the scheme of PDEs pdepe is able to solve.

pdepe is a solver for parabolic-elliptic PDEs that have the form described in the documentation (i.e. all equations have a second-order derivative term of the solution variable).

This is not true for both your first and your third equation.

You will have to discretize your equations in z-direction and solve the resulting system of ordinary differential equations using ODE15S, e.g.. Look up "method of lines" for more details.

[University Glasgow]

Okay, thank you.

[University Glasgow]

Please any idea on how to discretise mixed pdes?

[Torsten]

Do it the way you wanted to use pdepe: Define w = dtheta/dz. Then you have no mixed derivatives.

[University Glasgow]

Oh, okay, thank you.

[Torsten]

I think this system is not that easy to solve.

Where do the equations stem from ?

Are there discretization schemes in the literature especially designed for this system ?

[University Glasgow]

This stem from the Ericksen-Leslies equations linearised around the homeotropic alignment. I have been checking to see if I will get any discretization schemes in the literature but not yet found one.

[University Glasgow]

Below is the discretization that I tried and used 0de15s, but i recieved an error:

File: researhcase1.m Line: 73 Column: 58

Invalid expression. When calling a function or indexing a variable, use parentheses. Otherwise, check formismatched delimiters.

clear
clc
% parameters.
k1 = 6*10^(-12);    
eta1 = 0.0240;
apha3 = -0.001104;
gama1 = 0.1093;
d = 0.0002

d = 2.0000e-04

thetab = 0.0001

thetab = 1.0000e-04

% activity parameter
xi = 0.1; 
%TIME SPAN
Ts = 60.0;
tt = 0:100:2400;
t = tt/Ts;
%STEP LENGTH
L = 1.0;
Nz = 20;
zz= linspace(0, L, Nz);
z = zz;
dz = z(2) -z(1);
%IC
IC_A = zeros(1, Nz);
IC_B = zeros(1, Nz);
IC = [IC_A, IC_B];
%Solve the pde using ode15s
[t, u] = ode15s(@(t,y)f(t, y, u, v, Nz, k1, eta1, apha3, gama1, d, thetab, xi, dz), t, IC);

Unrecognized function or variable 'u'.

Error in solution (line 27)
[t, u] = ode15s(@(t,y)f(t, y, u, v, Nz, k1, eta1, apha3, gama1, d, thetab, xi, dz), t, IC);

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);

% Extract solutions 
U = u(:, 1:Nz);
V = v(:, 1:Nz);
W = v(:, Nz+1 : 2*Nz); 
% Reinput BC
U(:, 1) = 0;
V(:, 1) = 0; 
U(:, end) = (4*U(:, Nz-1)-U(:, Nz-2))./3;
%plot
figure(1)
plot(t, U(:, end), '--')

function duvdt = f(t, u, v, Nz, k1, eta1, apha3, gama1, d, thetab, xi, dz)
duvdt = zeros(length(u), 1);
dudt = zeros(Nz, 1);
dwdt = zeros(Nz, 1);
%define value
U = u(1:Nz);
V = v(1:Nz);
W = v(Nz+1 : 2*Nz); 
%BC
U(1) = 0;
V(1) = 0; 
U(end) = (4*U(Nz-1)-U(Nz-2))./3;
%interior for-loop 
for i=2:Nz-1
    w(i) = (U(i+1) - U(i-1))./2./dz;
    dvdz(i) = (V(i+1) - V(i-1))./2./dz;
    dv2dz(i) = (V(i+1) - 2*V(i) + V(i-1))/(dz*dz);
    du2dz(i) = (U(i+1) - 2*U(i) + U(i-1))/(dz*dz);
    dudt(i) = (1./gama1).*(k1.*du2dz(i) + apha3.*dvdz(i));
    dwdt(i) = (-1./apha3).*(eta1.* dv2dz(i) -w(i).*xi);
end
duvdz = [dudt; dwdt];
end

[Torsten]

Check the values for L and Ts.

The solution variables must be stored in one big vector (it's called y in the code above). You cannot supply u and v separately. So remove the parameters u and v from the list of input parameters to f. It seems to me that your vector of initial conditions has 2*Nz elements, but that you work with 3*Nz solution variables in f. Consequently, you also have to supply 3*Nz initial values and return duvdz as a vector of 3*Nz elements, not 2*Nz.

[University Glasgow]

I did what you said but still having the same problem. Please help out, I'm new to Matlab:

clear

clc

% parameters.

k1 = 6*10^(-12);

eta1 = 0.0240;

Ts = 60;

apha3 = -0.001104;

gama1 = 0.1093;

d = 0.0002

thetab = 0.0001

% activity parameter

xi = 0.1;

%TIME SPAN

tt = 0:10:200;

t= tt/Ts;

%STEP LENGTH

Nz = 20;

z= linspace(0, d, Nz);

dz = z(2) -z(1);

%IC

IC_A = zeros(1, Nz);

IC_B = zeros(1, Nz);

IC = [IC_A, IC_B];

%Solve the pde using ode15s

[t, u] = ode15s(@(t,y)f(t, y, Nz, k1, eta1, apha3, gama1, d, thetab, xi, dz), t, IC);

% Extract solutions

U = u(:, 1:Nz);

V = v(:, 1:Nz);

W = v(:, Nz+1 : 3*Nz);

% Reinput BC

U(:, 1) = 0;

V(:, 1) = 0;

U(:, end) = (4*U(:, Nz-1)-U(:, Nz-2))./3;

%plot

figure(1)

plot(t, U(:, end), '--')

function duvdt = f(t, Nz, k1, eta1, apha3, gama1, d, thetab, xi, dz)

duvdt = zeros(length(u), 1);

dudt = zeros(Nz, 1);

dwdt = zeros(Nz, 1);

%define value

U = u(1:Nz);

V = v(1:Nz);

W = v(Nz+1 : 3*Nz);

%BC

U(1) = 0;

V(1) = 0;

U(end) = (4*U(Nz-1)-U(Nz-2))./3;

%interior for-loop

for i=2:Nz-1

w(i) = (U(i+1) - U(i-1))./2./dz;

dvdz(i) = (V(i+1) - V(i-1))./2./dz;

dv2dz(i) = (V(i+1) - 2*V(i) + V(i-1))/dz*dz;

du2dz(i) = (U(i+1) - 2*U(i) + U(i-1))./dz*dz;

dudt(i) = (1./gama1).*(k1.*du2dz(i) + apha3.*dvdz(i));

dwdt(i) = (-1./apha3).*(eta1.* dv2dz(i) -w(i).*xi);

end

duvdz = [dudt; dwdt];

end

[Torsten]

The code below will not work, but it will at least show you what to do, I guess.

And don't repeat the error in the calculation of the second derivatives - I already corrected it above in your code and did it again below:

dv2dz(i) = (V(i+1) - 2*V(i) + V(i-1))/(dz*dz);
du2dz(i) = (U(i+1) - 2*U(i) + U(i-1))/(dz*dz);

instead of

dv2dz(i) = (V(i+1) - 2*V(i) + V(i-1))/dz*dz;
du2dz(i) = (U(i+1) - 2*U(i) + U(i-1))./dz*dz;

clear
clc
% parameters.
k1 = 6*10^(-12);    
eta1 = 0.0240;
Ts = 60;
apha3 = -0.001104;
gama1 = 0.1093;
d = 0.0002
thetab = 0.0001
% activity parameter
xi = 0.1; 
%TIME SPAN
tt = 0:10:200;
t= tt/Ts;
%STEP LENGTH
Nz = 20;
z= linspace(0, d, Nz);
dz = z(2) -z(1);
%IC
IC_A = zeros(1, Nz);    % Initial conditions for U
IC_B = zeros(1, Nz);    % Initial conditions for V
IC_C = zeros(1, Nz);    % Initial conditions for W
IC = [IC_A, IC_B, IC_C];
%Solve the pde using ode15s
[T,Y] = ode15s(@(t,y)f(t, y, Nz, k1, eta1, apha3, gama1, d, thetab, xi, dz), t, IC);
% Extract solutions 
U = Y(:, 1 : Nz);
V = Y(:, Nz+1 : 2*Nz);
W = Y(:, 2*Nz+1 : 3*Nz); 
% Reinput BC
U(:, 1) = 0;
V(:, 1) = 0; 
U(:, end) = (4*U(:, Nz-1)-U(:, Nz-2))./3;
%plot
figure(1)
plot(t, U(:, end), '--')
function duvwdt = f(t, y, Nz, k1, eta1, apha3, gama1, d, thetab, xi, dz)
duvwdt = zeros(length(y), 1);
dudt = zeros(Nz, 1);
dvdt = zeros(Nz ,1);
dwdt = zeros(Nz, 1);
%define value
U = y(1:Nz);
V = y(Nz+1:2*Nz);
W = y(2*Nz+1 : 3*Nz); 
%BC
U(1) = 0;
V(1) = 0; 
U(end) = (4*U(Nz-1)-U(Nz-2))./3;
%interior for-loop 
for i=2:Nz-1
    w(i) = (U(i+1) - U(i-1))./2./dz;
    dvdz(i) = (V(i+1) - V(i-1))./2./dz;
    dv2dz(i) = (V(i+1) - 2*V(i) + V(i-1))/(dz*dz);
    du2dz(i) = (U(i+1) - 2*U(i) + U(i-1))/(dz*dz);
    dudt(i) = (1./gama1).*(k1.*du2dz(i) + apha3.*dvdz(i));
    dvdt(i) = ...?
        dwdt(i) = (-1./apha3).*(eta1.* dv2dz(i) -w(i).*xi);
end
duvwdt = [dudt; dvdt ; dwdt];
end

[University Glasgow]

Okay, thank you

[Torsten]

I don't know if it's reasonable to solve the equations the way you try. I only changed the code that it might formally work.

You should search the literature for numerical approaches to solve the equations - standard approaches like yours will fail, I fear.

[University Glasgow]

Sure I will, thank you.