Please help me solve this material balance

[Sam Chak]

Hi @Rachael

Is it common to solve the material balance equations with a number of unknown parameters in your line of work? Or do you want to find the analytical solution for this Boundary Value Problem?

[Rachael]

[Rachael]

@Sam Chak If this would help solving the balance

[Sam Chak]

@Rachael,

You can follow the example of solving a system of PDEs in this link:

https://www.mathworks.com/help/matlab/math/solve-system-of-pdes.html

Use the template below to define the parameters and code the equations in the 'Local functions' section. If there is an error message, show the code.

x   = [0 0.005 0.01 0.05 0.1 0.2 0.5 0.7 0.9 0.95 0.99 0.995 1];
t   = [0 0.005 0.01 0.05 0.1 0.5 1 1.5 2];

m   = 0;
sol = pdepe(m, @pdefun, @pdeic, @pdebc, x, t);

u1  = sol(:,:,1);
u2  = sol(:,:,2);

surf(x,t,u1)
title('u_1(x,t)')
xlabel('Distance x')
ylabel('Time t')

surf(x,t,u2)
title('u_2(x,t)')
xlabel('Distance x')
ylabel('Time t')

%% -------------- Local functions -------------- 
% Partial Differential Equation to solve
function [c, f, s] = pdefun(x, t, u, dudx)
    c = [1; 1];
    f = [0.024; 0.17].*dudx;
    y = u(1) - u(2);
    F = exp(5.73*y) - exp(- 11.47*y);
    s = [-F; F];
end

% Initial Conditions
function u0 = pdeic(x) 
    u0 = [1; 0];
end

% Boundary Conditions
function [pl, ql, pr, qr] = pdebc(xl, ul, xr, ur, t)
    pl = [0; ul(2)];
    ql = [1; 0];
    pr = [ur(1)-1; 0];
    qr = [0; 1];
end

[Rachael]

Unrecognized function or variable 'pdeic'.

Error in pdepe (line 229)

temp = feval(ic,xmesh(1),varargin{:});

[Sam Chak]

@Rachael, click this icon and then copy/paste the entire MATLAB code to the grey field. Click this icon and see if it generate the "Unrecognized function" error message. The "pdepe" method is suggested by @Torsten.

[Rachael]

x = [18.4 15.7 22.6 21.8 29.8 21.1 21.9 17.9 21 16.7 16.2 22.1 19.1];

t = [0 24 48 72 96 120 144 168 192 216 240 264 288];

m = 0;

sol = pdepe(m, @pdefun, @pdeic, @pdebc, x, t);

u1 = sol(:,:,1);

u2 = sol(:,:,2);

surf(x,t,u1)

title('u_1(x,t)')

xlabel('Distance x')

ylabel('Time t')

surf(x,t,u2)

title('u_2(x,t)')

xlabel('Distance x')

ylabel('Time t')

%% -------------- Local functions --------------

% Partial Differential Equation to solve

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

c = [1; 1];

f = [0.00000256; 0.0000198].*dudx;

y = u(1) - u(2);

F = 0.0017*exp(-0.082*y) - -5.63*exp(-0.082*y);

s = [-F; F];

end

% Initial Conditions

function u0 = pdeic(x)

u0 = [1; 0];

end

% Boundary Conditions

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

pl = [0; ul(2)];

ql = [1; 0];

pr = [ur(1)-1; 0];

qr = [0; 1];

end

Error using pdepe

The entries of XMESH must be strictly increasing.

[Rachael]

x   = [0.05 0.1 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.61];
t   = [0 24 48 72 96 120 144 168 192 216 240 264 288];

m   = 0;
sol = pdepe(m, @pdefun, @pdeic, @pdebc, x, t);

u1  = sol(:,:,1);
u2  = sol(:,:,2);

surf(x,t,u1)
title('u_1(x,t)')
xlabel('Distance x')
ylabel('Time t')

surf(x,t,u2)
title('u_2(x,t)')
xlabel('Distance x')
ylabel('Time t')

%% -------------- Local functions -------------- 
% Partial Differential Equation to solve

function [c, f, s] = pdefun(x, t, u, dudx)
    c = [1; 1];
    f = [0.00000256; 0.0000198].*dudx;
    y = u(1) - u(2);
    F = 0.0017*exp(-0.082*y) - -5.63*exp(-0.082*y);
    s = [-F; F];
end

% Initial Conditions

function u0 = pdeic(x) 
    u0 = [1; 0];
end

% Boundary Conditions

function [pl, ql, pr, qr] = pdebc(xl, ul, xr, ur, t)
    pl = [0; ul(2)];
    ql = [1; 0];
    pr = [ur(1)-1; 0];
    qr = [0; 1];
end
Unrecognized function or variable 'pdeic'.

Function definitions in a script must appear at the end of the file.
Move all statements after the "pdebc" function definition to before the first local function definition.

Error in pdepe (line 229)
temp = feval(ic,xmesh(1),varargin{:});

Answer 1

[Torsten]

This is a typical problem for MATLAB's "pdepe". So I suggest you use this integrator for partial differential equations.

[Rachael]

x = [18.4 15.7 22.6 21.8 29.8 21.1 21.9 17.9 21 16.7 16.2 22.1 19.1];

t = [0 24 48 72 96 120 144 168 192 216 240 264 288];

m = 0;

sol = pdepe(m, @pdefun, @pdeic, @pdebc, x, t);

u1 = sol(:,:,1);

u2 = sol(:,:,2);

surf(x,t,u1)

title('u_1(x,t)')

xlabel('Distance x')

ylabel('Time t')

surf(x,t,u2)

title('u_2(x,t)')

xlabel('Distance x')

ylabel('Time t')

%% -------------- Local functions --------------

% Partial Differential Equation to solve

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

c = [1; 1];

f = [0.00000256; 0.0000198].*dudx;

y = u(1) - u(2);

F = 0.0017*exp(-0.082*y) - -5.63*exp(-0.082*y);

s = [-F; F];

end

% Initial Conditions

function u0 = pdeic(x)

u0 = [1; 0];

end

% Boundary Conditions

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

pl = [0; ul(2)];

ql = [1; 0];

pr = [ur(1)-1; 0];

qr = [0; 1];

end

Error using pdepe

The entries of XMESH must be strictly increasing.

[Torsten]

x = [18.4 15.7 22.6 21.8 29.8 21.1 21.9 17.9 21 16.7 16.2 22.1 19.1];

This should be the spatial discretization points of your x-interval of integration (most probably of the column height). Do you really want to integrate from x = 15.7 to x = 29.8 ? Then input this vector in increasing order to "pdepe".

[Rachael]

x = [0.05 0.1 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.61];

t = [0 24 48 72 96 120 144 168 192 216 240 264 288];

m = 0;

sol = pdepe(m, @pdefun, @pdeic, @pdebc, x, t);

u1 = sol(:,:,1);

u2 = sol(:,:,2);

surf(x,t,u1)

title('u_1(x,t)')

xlabel('Distance x')

ylabel('Time t')

surf(x,t,u2)

title('u_2(x,t)')

xlabel('Distance x')

ylabel('Time t')

%% -------------- Local functions --------------

% Partial Differential Equation to solve

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

c = [1; 1];

f = [0.00000256; 0.0000198].*dudx;

y = u(1) - u(2);

F = 0.0017*exp(-0.082*y) - -5.63*exp(-0.082*y);

s = [-F; F];

end

% Initial Conditions

function u0 = pdeic(x)

u0 = [1; 0];

end

% Boundary Conditions

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

pl = [0; ul(2)];

ql = [1; 0];

pr = [ur(1)-1; 0];

qr = [0; 1];

end

Unrecognized function or variable 'pdeic'.

Error in pdepe (line 229)

temp = feval(ic,xmesh(1),varargin{:});

[Rachael]

Its still showing some error

[Torsten]

@Sam Chak supplied a working code from the pdepe examples collection. You made changes to the code - so there is no guarantee that the code will work thereafter. In your case, there are problems for t > 1.17. So reduce the time span and see what happens short before the solver quits. Obviously, the solution blows up.

x   = [0.05 0.1 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.61];
t   = [0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.17];
m   = 0;
sol = pdepe(m, @pdefun, @pdeic, @pdebc, x, t);
u1  = sol(:,:,1);
u2  = sol(:,:,2);
surf(x,t,u1)
title('u_1(x,t)')
xlabel('Distance x')
ylabel('Time t')

surf(x,t,u2)
title('u_2(x,t)')
xlabel('Distance x')
ylabel('Time t')

%% -------------- Local functions -------------- 
% Partial Differential Equation to solve
function [c, f, s] = pdefun(x, t, u, dudx)
    c = [1; 1];
    f = [0.00000256; 0.0000198].*dudx;
    y = u(1) - u(2);
    F = 0.0017*exp(-0.082*y) - -5.63*exp(-0.082*y);
    s = [-F; F];
end
% Initial Conditions
function u0 = pdeic(x) 
    u0 = [1; 0];
end
% Boundary Conditions
function [pl, ql, pr, qr] = pdebc(xl, ul, xr, ur, t)
    pl = [0; ul(2)];
    ql = [1; 0];
    pr = [ur(1)-1; 0];
    qr = [0; 1];
end

[Rachael]

x   = [0.05 0.1 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.61];
t   = [0 24 48 72 96 120 144 168 192 216 240 264 288];

m   = 0;
sol = pdepe(m, @pdefun, @icfun, @bcfun, x, t);

Error using solution>pdefun
Too many input arguments.

Error in pdepe (line 246)
[c,f,s] = feval(pde,xi(1),t(1),U,Ux,varargin{:});

u1  = sol(:,:,1);
u2  = sol(:,:,2);

surf(x,t,u1)
title('u_1(x,t)')
xlabel('Distance x')
ylabel('Time t')

surf(x,t,u2)
title('u_2(x,t)')
xlabel('Distance x')
ylabel('Time t')

%% -------------- Local functions -------------- 
% Partial Differential Equation to solve

function [c, f, s] = pdefun(x, t, dudx)
D1= 0.00000190;
D2=0.0000198;
p=1;
v1=0.0011;
v2=0.0026;
r=0.000006;
k=0.082;
C=0.0205;
R=8.205;
T=305.15;
M=16;
e=0.036;

c = [1; 1];
f = [D1*dudx-v1; p*D2*dudx-v2*p].*dudx;
F = -r+k*(C-c(x,t));
J =(-(R*T*e)/M)+k*(C-c(x,t));
s = [F; J];
end

% Initial Conditions

function u0 = icfun(x) 
u0 = [x; 0];  
end

% Boundary Conditions

function [pl, ql, pr, qr] = bcfun(xl, ul, xr, ur, t)
D1= 0.00000190;
D2=0.0000198;
v1=0.0011;
v2=0.0026;

pl = [0; (v1/D1)*c(0,t)];
ql = [0; 0];
pr = [(v2/D2)*c(0,t); 0];
qr = [0; 0];
end

[Torsten]

Please use the same variable names as in your equation sheet.

The definition of f and your boundary condition settings are definitely wrong.

And pdefun has 4 inputs, not 3.

[Rachael]

x   = [0.05 0.1 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.61];
t   = [0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.17];

m   = 0;
sol = pdepe(m, @pdefun, @icfun, @bcfun, x, t);

Error using pdepe
Unexpected output of PDEFUN. For this problem PDEFUN must return three column vectors of length 2.

u1  = sol(:,:,1);
u2  = sol(:,:,2);

surf(x,t,u1)
title('u_1(x,t)')
xlabel('Distance x')
ylabel('Time t')

surf(x,t,u2)
title('u_2(x,t)')
xlabel('Distance x')
ylabel('Time t')

%% -------------- Local functions -------------- 
% Partial Differential Equation to solve

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

D1= 0.00000190;
D2=0.0000198;
p=1;
v1=0.0011;
v2=0.0026;
r=0.000006;
k=0.082;
C=0.0205;
R=8.205;
T=305.15;
M=16;
e=0.036;

c = [1; 1];
f = [D1-v1; p*(D2-v2)].*dudx;
F = -r+k*(C-c);
J =(-(R*T*e)/M)+k*(C-c);
s = [F; J];
end

% Initial Conditions

function u0 = icfun(x) 
u0 = [x; 0];  
end

% Boundary Conditions

function [pl, ql, pr, qr] = bcfun(xl, ul, xr, ur, t)
pl = [0; ul(0.61)];
ql = [0; 0];
pr = [0; ur(0.61)];
qr = [0; 0];
end

[Torsten]

x   = [0.05 0.1 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.61];
t   = [0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.17];

m   = 0;
sol = pdepe(m, @pdefun, @icfun, @bcfun, x, t);

u1  = sol(:,:,1);
u2  = sol(:,:,2);

surf(x,t,u1)
title('u_1(x,t)')
xlabel('Distance x')
ylabel('Time t')

surf(x,t,u2)
title('u_2(x,t)')
xlabel('Distance x')
ylabel('Time t')

%% -------------- Local functions -------------- 
% Partial Differential Equation to solve

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

D1= 0.00000190;
D2=0.0000198;
p=1;
v1=0.0011;
v2=0.0026;
r=0.000006;
k=0.082;
C=0.0205;
R=8.205;
T=305.15;
M=16;
e=0.036;

c = [1; 1];
f = [D1; p*D2].*dudx;
F = -r+k*(C-u(1));
J =(-(R*T*e)/M)+k*(C-u(1));
s = [F-v1*dudx(1); J-p*v2*dudx(2)];
end

% Initial Conditions

function u0 = icfun(x) 
u0 = [x; 0];  
end

% Boundary Conditions

function [pl, ql, pr, qr] = bcfun(xl, ul, xr, ur, t)
v1=0.0011;
v2=0.0026;
C=0.0205;
pl = [-v1*(ul(1)-C); -v2*(ul(2)-C)];
ql = [1; 1];
pr = [0; 0];
qr = [1; 1];
end

Answer 2

[Sam Chak]

Hi @Rachael

@Torsten has shown you how to solve PDEs in MATLAB. From the PDEs,

you need to rearrange and rewrite the equations in the form that the pdepe() solver expects: c, f, s.

Only then, you can code the local function pdefun(x, t, u, dudx).