System of differential equation. How to solve?

Hi to all,

it's my first time I'm trying to solve differential equations with MATLAB (2018a). I'm trying to solve a three differential equations system as the following one

%Equations
eqElec = Ce(t,z)*diff(Te(t,z),t)-diff(ke(t,z)*diff(Te(t,z),z),z) == -Gl(t,z)*(Te(t,z)-Ti(t,z)) + sL(t,z); 
eqLat  = Ci*diff(Ti(t,z),t)-diff(ki*diff(Ti(t,z),z),z) == Gl(t,z)*(Te(t,z)-Ti(t,z));
eqDiff = diff(ne(t,z),t) == DD(t,z)*diff(ne(t,z),z,2);

It involves two systems that exchange heat. The boundary conditions are

%Boundary conditions
cond1 = Te(t,zf) == T0;
cond2 = Te(t0,z) == T0;
cond3 = Ti(t,zf) == T0;
cond4 = Ti(t0,z) == T0;
cond5 = ne(t0,z) == ne0;
cond6 = ne(t,zf) == ne0;
cond7 = ki*diff(Ti(t,z0),z) == 0;
cond8 = ke(t,z0)*diff(Te(t,z0),z) == 0;
cond9 = DD(t,z0)*diff(ne(t,z0),z) == 0;

All the main functions have been declared as symbolic objects

%Declare symbolic functions
syms Te(t,z) Ti(t,z) ne(t,z) Ve(t,z) aa(t,z) sL(t,z) Ce(t,z) fei(t,z) ke(t,z) se(t,z) ue(t,z) DD(t,z) Gl(t,z) lfp(t,z) Jth(t,z) Jph(t,z);

How can I solve this system? Probably it can be resolved only numerically but I don't know whick solver I have to use. ode45?

Thanks in advance!

Answer 1

[Amish]

Hi Riccardo,

In order to solve a system of partial differential equations (PDEs) in MATLAB, especially with boundary conditions, you typically need to use numerical solvers designed for PDEs. In MATLAB, pdepe is a function that can handle systems of parabolic and elliptic PDEs in one spatial dimension and time.

A general approach to solve PDE involves, converting the equations into the form expected by pdepe function, followed by specifying the initial and boundary conditions in the correct format and then using the function to solve the equations. We can then plot the solutions after extracting them.

Below is a generic boiler plate code for solving PDEs:

% Define the spatial domain
z0 = 0; % Start of spatial domain
zf = 1; % End of spatial domain
m = 0;  % Symmetry parameter (0 for slab, 1 for cylindrical, 2 for spherical)

% Define the time domain and time points
t0 = 0;
tf = 10; 
tspan = linspace(t0, tf, 100); 

% Define the PDE system
function [c, f, s] = pdefun(z, t, u, DuDz)
    % u(1) = Te, u(2) = Ti, u(3) = ne
    % DuDz(1) = dTe/dz, DuDz(2) = dTi/dz, DuDz(3) = dne/dz

    % Coefficients for the PDEs
    Ce_val = 1; 
    ke_val = 1;
    Gl_val = 1;
    sL_val = 0;
    Ci_val = 1;
    ki_val = 1;
    DD_val = 1;

    % Define the system of equations
    c = [Ce_val; Ci_val; 1];
    f = [ke_val * DuDz(1); ki_val * DuDz(2); DD_val * DuDz(3)];
    s = [-Gl_val * (u(1) - u(2)) + sL_val; Gl_val * (u(1) - u(2)); 0];
end

% Initial conditions
function u0 = icfun(z)
    T0 = 300; 
    ne0 = 1;  
    u0 = [T0; T0; ne0];
end

% Boundary conditions
function [pl, ql, pr, qr] = bcfun(zl, ul, zr, ur, t)
    T0 = 300; 
    ne0 = 1; 

    % Left boundary conditions (z = z0)
    pl = [0; 0; 0];
    ql = [1; 1; 1];

    % Right boundary conditions (z = zf)
    pr = [ul(1) - T0; ur(2) - T0; ur(3) - ne0];
    qr = [0; 0; 0];
end

% Solve
sol = pdepe(m, @pdefun, @icfun, @bcfun, linspace(z0, zf, 100), tspan);

% Extract solutions
Te_sol = sol(:,:,1);
Ti_sol = sol(:,:,2);
ne_sol = sol(:,:,3);

% Plot the results
figure;
subplot(3,1,1);
surf(linspace(z0, zf, 100), tspan, Te_sol);
title('Te');
xlabel('z');
ylabel('t');
zlabel('Te');

subplot(3,1,2);
surf(linspace(z0, zf, 100), tspan, Ti_sol);
title('Ti');
xlabel('z');
ylabel('t');
zlabel('Ti');

subplot(3,1,3);
surf(linspace(z0, zf, 100), tspan, ne_sol);
title('ne');
xlabel('z');
ylabel('t');
zlabel('ne');

For more information and additional examples, you can visit the following documentation:

pdepe: https://www.mathworks.com/help/matlab/ref/pdepe.html

The following community post might also be helpful:

https://www.mathworks.com/matlabcentral/answers/2081036-how-to-output-quantities-within-pdefun-bcfun-icfun-in-pdepe

Hope this helps!