Why am I getting error using assignin and syms while solving pde using pdepe

So I am trying to solve the non dimensionalized form of an equation which describes pore diffusion resistance combined with surface kinetics. So the equation is a partial differential equation and its boundary condition is a system of coupled ODEs. I have attached the file regarding the equation above. Refer to the non dimensional section for the equations. Here is the code that I have written which uses pdepe and dsolve.

function coupled_pde_ode()
    % Define parameters
    epsilon_p = input('Porosity (epsilon_p): ');
    De = input('Effective diffusivity (De, in m^2/s): ');
    R = input('Radius of the pellet (R, in m): ');
    kf = input('Mass transfer coefficient (kf, in m/s): ');
    kR = input('Reaction rate constant (kR, in 1/s): ');
    C0 = input('Initial concentration in pellet pores (C0, arbitrary units): ');
    mp = input('Weight of catalyst beads (mp, in kg): ');
    rho_p = input('Overall density of catalyst pellets (rho_p, in kg/m^3): ');
    V = input('Volume of reactor (V, in m^3): ');
    
    % Derived non-dimensional parameters
    phi = (R / 3) * sqrt(kR / De); % Thiele modulus
    Bi = kf * R / De; % Biot number
    alpha = (3 * epsilon_p * mp) / (V * rho_p); % Coupling parameter
    tau_end = 10; % End time in dimensionless units

    % Non-dimensionalized spatial and temporal discretization
    rbar = linspace(0, 1, 50); % Non-dimensional radius (0 <= rbar <= 1)
    tau = linspace(0, tau_end, 100); % Non-dimensional time

    % Solve PDE for Cpbar using pdepe
    m = 2; % Spherical symmetry
    sol_pde = pdepe(m, ...
                    @(rbar, tau, Cpbar, dCpbar_drbar) pde_eqn(rbar, tau, Cpbar, dCpbar_drbar, epsilon_p, phi), ...
                    @initial_condition, ...
                    @(xl, CpbarL, xr, CpbarR, tau) boundary_conditions(xl, CpbarL, xr, CpbarR, tau, Bi, phi), ...
                    rbar, tau);
    
    % Extract Cpbar at rbar = 1
    Cpbar_surface = sol_pde(:, end); % Concentration at rbar = 1 over time
    
    % Solve coupled ODE for Cbbar using dsolve
    syms Cbbar(tau)
    dCbbar_dtau = -alpha * (Bi * (Cbbar - Cpbar_surface') - 3 * phi^2 * Cpbar_surface');
    Cbbar_sol = dsolve(dCbbar_dtau, Cbbar(0) == 1);
    
    % Evaluate Cbbar over time (not plotted)
    tau_values = tau; % Time grid
    Cbbar_values = double(subs(Cbbar_sol, tau, tau_values));

    % Plot results: Concentration vs. Radius
    figure;
    for i = 1:10:length(tau)
        plot(rbar, sol_pde(i, :), 'DisplayName', sprintf('\\tau = %.1f', tau(i)));
        hold on;
    end
    xlabel('\bar{r} (Dimensionless Radius)');
    ylabel('\bar{C}_p (Dimensionless Pore Concentration)');
    legend;
    title('Pore Concentration Profiles (\bar{C}_p) vs. Radius');
    grid on;
    hold off;

    % Nested functions
    function [c, f, s] = pde_eqn(rbar, tau, Cpbar, dCpbar_drbar, epsilon_p, phi)
        % Non-dimensional PDE coefficients
        c = epsilon_p; % Time-dependent term coefficient
        f = dCpbar_drbar; % Diffusion term coefficient
        s = -9 * phi^2 * Cpbar; % Source term (reaction rate)
    end

    function Cpbar0 = initial_condition(rbar)
        % Non-dimensional initial condition
        Cpbar0 = ones(size(rbar)); % Initial concentration equals C0 (normalized to 1)
    end

    function [pl, ql, pr, qr] = boundary_conditions(xl, CpbarL, xr, CpbarR, tau, Bi, phi)
        % Boundary conditions for the non-dimensional problem
        % At rbar = 0 (center of the sphere), Neumann BC: dCpbar/drbar = 0
        pl = 0;
        ql = 1;

        % At rbar = 1 (outer surface of the sphere), flux term
        pr = Bi * (1 - CpbarR) - 3 * phi^2 * CpbarR;
        qr = 1;
    end
end

Here are the values I entered

Porosity (epsilon_p):

0.5

Effective diffusivity (De, in m^2/s):

1.5e-10

Radius of the pellet (R, in m):

1e-3

Mass transfer coefficient (kf, in m/s):

2.8e-3

Reaction rate constant (kR, in 1/s):

5e-12

Initial concentration in pellet pores (C0, arbitrary units):

1

Weight of catalyst beads (mp, in kg):

0.1e-3

Overall density of catalyst pellets (rho_p, in kg/m^3):

850

Volume of reactor (V, in m^3):

0.1

Here is the error message I'm receiving

Error using assignin

Attempt to add "Cbbar" to a static workspace.

Error in syms (line 331)

assignin('caller', name, xsym);

^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

Error in Testing (line 35)

syms Cbbar(tau)

Answer 1

[Torsten]

Here is the code for the dimensional model.

It should be no problem to rewrite it in the non-dimensional form.

Note that the sign in S1c of your attached document is incorrect:

it should be

kf*(Cp-Cb)

instead of

kf*(Cb-Cp)

% Define parameters
epsilon_p = 0.5;%input('Porosity (epsilon_p): ');
De = 1.5e-10;%input('Effective diffusivity (De, in m^2/s): ');
R = 1e-3;%input('Radius of the pellet (R, in m): ');
kf = 2.8e-3;%input('Mass transfer coefficient (kf, in m/s): ');
kR = 5.0e-12;%input('Reaction rate constant (kR, in 1/s): ');
C0 = 1;%input('Initial concentration in pellet pores (C0, arbitrary units): ');
mp = 1e-4;%input('Weight of catalyst beads (mp, in kg): ');
rho_p = 850;%input('Overall density of catalyst pellets (rho_p, in kg/m^3): ');
V = 1e-1;%input('Volume of reactor (V, in m^3): ');
    
% spatial and temporal discretization
nr = 100;
r = linspace(0, R, nr).'; % radius (0 <= rbar <= 1)
dr = r(2) - r(1);
tend = 100; % End time 
nt = 10;
tspan = linspace(0, tend, nt); % time
    
% Initial conditions
Cp0 = zeros(nr,1);
Cb0 = 1.0;
y0 = [Cp0;Cb0];
    
[T,Y] = ode15s(@(t,y)fun(t,y,nr,r,dr,epsilon_p,De,R,kf,kR,C0,mp,rho_p,V),tspan,y0);
   
Cp = Y(:,1:end-1);
Cb = Y(:,end);

figure(1)
plot(r,Cp)

figure(2)
plot(T,Cb)

    
function dydt = fun(t,y,nr,r,dr,epsilon_p,De,R,kf,kR,C0,mp,rho_p,V)
  Cp = y(1:end-1);
  Cb = y(end);
  dCp_dt = zeros(nr,1);
  dCp_dt(1) = 3/((r(1)+dr/2)^3-r(1)^3)*De*...
              ((r(1)+dr/2)^2*(Cp(2)-Cp(1))/dr)-...
              kR * Cp(1);
  dCp_dt(2:nr-1) = 3./((r(2:nr-1)+dr/2).^3-(r(2:nr-1)-dr/2).^3)*De.*...
                   ((r(2:nr-1)+dr/2).^2.*(Cp(3:nr)-Cp(2:nr-1))/dr-...
                    (r(2:nr-1)-dr/2).^2.*(Cp(2:nr-1)-Cp(1:nr-2))/dr)-...
                   kR * Cp(2:nr-1);
  dCp_drR = -(kf*(Cp(nr)-Cb) - R/3*kR*Cp(nr))/De;
  dCp_dt(nr) = 3/(r(nr)^3-(r(nr)-dr/2)^3)*De*...
               (r(nr)^2*dCp_drR - ...
                (r(nr)-dr/2)^2*(Cp(nr)-Cp(nr-1))/dr)-...
                kR * Cp(nr);
  dCp_dt = dCp_dt/epsilon_p;
  dCb_dt = -1/V * mp/rho_p * 3/R * De * dCp_drR;
  dydt = [dCp_dt;dCb_dt];
end

[Aditya]

Thanks for your time and effort. I'll definitely look into this code. But when you sent in your initial response I tried to use a different method. Rather than substituting the equation S4 into S3c I tried a numerical method. I took a guess value of dCpbar/drbar at rbar=1. Calculated Cbbar from S4. Then solved S3 using Cbbar and dCpbar/drbar. Then I kept updating the values till they converge at abs(Cpbar(tau) - Cpbar(tau-1))<=epsilon where epsilon is a very small value. Here is the code.

function iterative_pde_ode_solver()
    % Input parameters
    disp('Enter the following parameters:');
    epsilon_p = input('Porosity (epsilon_p): ');
    De = input('Effective diffusivity (De, in m^2/s): ');
    Rp = input('Radius of the pellet (Rp, in m): ');
    kf = input('Mass transfer coefficient (kf, in m/s): ');
    kR = input('Reaction rate constant (kR, in 1/s): ');
    C0 = input('Initial bulk concentration (C0, arbitrary units): ');
    mp = input('Weight of catalyst beads (mp, in kg): ');
    rhop = input('Density of catalyst pellet (rhop, in kg/m^3): ');
    V = input('Volume of reaction medium (V, in m^3): ');
    epsilon = input('Convergence threshold (epsilon, e.g., 0.0001): ');
    
    % Derived non-dimensional parameters
    phi = (Rp / 3) * sqrt(kR / De); % Thiele modulus
    Bi = kf * Rp / De; % Biot number
    alpha = (3 * epsilon_p * mp) / (V * rhop); % Coupling parameter
    rbar = linspace(0, 1, 50); % Non-dimensional radius
    tau = linspace(0, 10, 100); % Non-dimensional time
    
    % Initial guess for dCpbar/drbar at rbar = 1
    dCpbar_drbar_guess = 1; % Arbitrary initial guess
    
    % Initialize variables
    Cpbar_prev = ones(length(rbar), 1); % Initial condition for Cpbar
    converged = false;
    iteration = 0;
    
    % Iterative solver
    while ~converged
        iteration = iteration + 1;
        fprintf('Iteration %d\n', iteration);
        
        % Solve ODE for Cbbar
        [tau_sol, Cbbar_sol] = ode45(@(tau, Cbbar) ...
            ode_rhs(tau, Cbbar, alpha, dCpbar_drbar_guess), tau, 1);
        
        % Solve PDE for Cpbar
        m = 2; % Spherical symmetry
        sol = pdepe(m, ...
                    @(rbar, tau, Cpbar, dCpbar_drbar) ...
                        pde_eqn(rbar, tau, Cpbar, dCpbar_drbar, epsilon_p, phi), ...
                    @initial_condition, ...
                    @(xl, CpbarL, xr, CpbarR, tau) ...
                        boundary_conditions(xl, CpbarL, xr, CpbarR, tau, Bi, phi, Cbbar_sol, tau_sol), ...
                    rbar, tau);
        
        % Extract new Cpbar
        Cpbar_new = sol(end, :)'; % Pore concentration at final time step
        
        % Convergence check
        if max(abs(Cpbar_new - Cpbar_prev)) <= epsilon
            converged = true;
        else
            % Update flux guess
            dCpbar_drbar_guess = (Cpbar_new(end) - Cpbar_new(end-1)) / (rbar(end) - rbar(end-1));
            Cpbar_prev = Cpbar_new; % Update previous Cpbar
        end
    end
    
    % Plot results
    figure;
    plot(rbar, Cpbar_new, 'LineWidth', 2);
    xlabel('\bar{r} (Dimensionless Radius)');
    ylabel('\bar{C}_p (Dimensionless Pore Concentration)');
    title('Final Pore Concentration Profile (\bar{C}_p) vs Radius');
    grid on;
    
    fprintf('Converged after %d iterations.\n', iteration);
    
    % Nested Functions
    function dCbbar_dt = ode_rhs(tau, Cbbar, alpha, dCpbar_drbar)
        dCbbar_dt = -alpha * dCpbar_drbar; % ODE for Cbbar
    end

    function [c, f, s] = pde_eqn(rbar, tau, Cpbar, dCpbar_drbar, epsilon_p, phi)
        c = epsilon_p; % Time-dependent term coefficient
        f = dCpbar_drbar; % Diffusion term coefficient
        s = -3 * phi^2 * Cpbar; % Source term (reaction rate)
    end

    function Cpbar0 = initial_condition(rbar)
        Cpbar0 = ones(size(rbar)); % Initial condition for Cpbar
    end

    function [pl, ql, pr, qr] = boundary_conditions(xl, CpbarL, xr, CpbarR, tau, Bi, phi, Cbbar_sol, tau_sol)
        % Interpolate Cbbar to match the current time tau
        Cbbar = interp1(tau_sol, Cbbar_sol, tau, 'linear', 'extrap');
        
        % At rbar = 0 (symmetry condition)
        pl = 0;
        ql = 1;
        
        % At rbar = 1 (outer surface of the pellet)
        pr = Bi * (Cbbar - CpbarR) - 3 * phi^2 * CpbarR;
        qr = 1;
    end
end

I just wanted to know your thoughts on it. Although I'll still have to debugg the code since I'm getting the time integration failure error here too

Enter the following parameters:

Porosity (epsilon_p):

0.5

Effective diffusivity (De, in m^2/s):

1.5e-10

Radius of the pellet (Rp, in m):

1e-3

Mass transfer coefficient (kf, in m/s):

2.8e-3

Reaction rate constant (kR, in 1/s):

5e-2

Initial bulk concentration (C0, arbitrary units):

10

Weight of catalyst beads (mp, in kg):

0.1e-3

Density of catalyst pellet (rhop, in kg/m^3):

850

Volume of reaction medium (V, in m^3):

0.1

Convergence threshold (epsilon, e.g., 0.0001):

0.0001

Iteration 1

Warning: Failure at t=1.884177e-04. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (4.336809e-19) at time t.

> In ode15s (line 625)

In pdepe (line 290)

In Testing (line 41)

Warning: Time integration has failed. Solution is available at requested time points up to t=0.000000e+00.

> In pdepe (line 306)

In Testing (line 41)

Converged after 1 iterations.

Warning: Error in state of SceneNode.

String scalar or character vector must have valid interpreter syntax:

Final Pore Concentration Profile (\\bar{C}_p) vs Radius

> In defaulterrorcallback (line 12)

In matlab.ui.internal.FigureServices.getUniqueChannelId (line 106)

In mls.internal.figureCreated

In mls.internal.figureCreated

Warning: Error in state of SceneNode.

String scalar or character vector must have valid interpreter syntax:

\\bar{r} (Dimensionless Radius)

> In defaulterrorcallback (line 12)

In matlab.ui.internal.FigureServices.getUniqueChannelId (line 106)

In mls.internal.figureCreated

In mls.internal.figureCreated

Warning: Error in state of SceneNode.

String scalar or character vector must have valid interpreter syntax:

\\bar{C}_p (Dimensionless Pore Concentration)

> In defaulterrorcallback (line 12)

In matlab.ui.internal.FigureServices.getUniqueChannelId (line 106)

In mls.internal.figureCreated

In mls.internal.figureCreated

Once again I would really like to thank you for the help and support you are providing. I'll take a look into the method of lines method too.

[Torsten]

In order to solve (S4) over time, you should insert (S3c) into (S4). (As mentionned above, there is a sign error in (S3c) in the convective flux term.)

Then, from the solution of "pdepe", you have to extract the values for Cpbar at rbar = 1 over time and use interpolation in ode_rhs to the times "tau" requested by ode45 (you shouldn't name them "tau" in the list of input parameters as they differ from the tau-vector used as output times for "pdepe").

So all in all: you might succeed, but iteration usually is the last refuge. The cleanest way is to define the PDE and ODE together and use an ODE integrator that solves all the resulting equations simultaneously.

I think that

dCpbar_drbar_guess = (Cpbar_new(end) - Cpbar_new(end-1)) / (rbar(end) - rbar(end-1));

is incorrect since

Cpbar_new(end) - Cpbar_new(end-1)

looks like a difference in time, divided by

rbar(end) - rbar(end-1)

which is a difference in space.

[Walter Roberson]

You need

    xlabel('\bar{r} (Dimensionless Radius)', 'interpreter', 'latex');
    ylabel('\bar{C}_p (Dimensionless Pore Concentration)', 'interpreter', 'latex');
    title('Final Pore Concentration Profile (\bar{C}_p) vs Radius', 'interpreter', 'latex');

[Aditya]

Also suppose I had to update this code to consider nth order reaction kinetics (using method of lines). So would making these changes be enough.

% Define parameters
epsilon_p = 0.5; % Porosity
De = 1.5e-12; % Effective diffusivity (m^2/s)
R = 1e-3; % Radius of the pellet (m)
kf = 2.8e-3; % Mass transfer coefficient (m/s)
kR = 5.0e-12; % Reaction rate constant (1/s)
C0 = 1; % Initial concentration in the pellet pores
mp = 1e-4; % Weight of catalyst beads (kg)
rho_p = 850; % Density of catalyst pellets (kg/m^3)
V = 1e-1; % Volume of the reactor (m^3)
n = 1; % Order of reaction

% Spatial and temporal discretization
nr = 100; % Number of radial nodes
r = linspace(0, R, nr).'; % Radial positions
dr = r(2) - r(1); % Spatial step size
tend = 100; % End time (s)
nt = 10; % Number of time points
tspan = linspace(0, tend, nt); % Time discretization

% Initial conditions
Cp0 = zeros(nr, 1); % Initial pore concentration
Cb0 = 1.0; % Initial bulk concentration
y0 = [Cp0; Cb0]; % Combined state vector

% Solve the system using ode15s
[T, Y] = ode15s(@(t, y) fun(t, y, nr, r, dr, epsilon_p, De, R, kf, kR, C0, mp, rho_p, V, n), tspan, y0);

% Extract pore and bulk concentrations
Cp = Y(:, 1:end-1); % Pore concentration
Cb = Y(:, end); % Bulk concentration

% Plot results
figure(1);
plot(r, Cp);
title('Pore Concentration Profile over Radius');
xlabel('Radius (m)');
ylabel('Concentration in Pellet (C_p)');

figure(2);
plot(T, Cb);
title('Bulk Concentration vs Time');
xlabel('Time (s)');
ylabel('Bulk Concentration (C_b)');

% Function to calculate derivatives
function dydt = fun(t, y, nr, r, dr, epsilon_p, De, R, kf, kR, C0, mp, rho_p, V, n)
    Cp = y(1:end-1); % Pore concentration
    Cb = y(end); % Bulk concentration
    dCp_dt = zeros(nr, 1); % Initialize dCp_dt

    % At r = 0 (center of the sphere)
    dCp_dt(1) = 3 / ((r(1) + dr/2)^3 - r(1)^3) * De * ...
                ((r(1) + dr/2)^2 * (Cp(2) - Cp(1)) / dr) - ...
                kR * Cp(1)^n;

    % For interior nodes (0 < r < R)
    dCp_dt(2:nr-1) = 3 ./ ((r(2:nr-1) + dr/2).^3 - (r(2:nr-1) - dr/2).^3) .* De .* ...
                     ((r(2:nr-1) + dr/2).^2 .* (Cp(3:nr) - Cp(2:nr-1)) / dr - ...
                      (r(2:nr-1) - dr/2).^2 .* (Cp(2:nr-1) - Cp(1:nr-2)) / dr) - ...
                     kR * Cp(2:nr-1).^n;

    % At r = R (surface of the pellet)
    dCp_drR = -(kf * (Cp(nr) - Cb) - R/3 * kR * Cp(nr)^n) / De;
    dCp_dt(nr) = 3 / (r(nr)^3 - (r(nr) - dr/2)^3) * De * ...
                 (r(nr)^2 * dCp_drR - (r(nr) - dr/2)^2 * (Cp(nr) - Cp(nr-1)) / dr) - ...
                 kR * Cp(nr)^n;

    % Normalize by porosity
    dCp_dt = dCp_dt / epsilon_p;

    % Bulk concentration ODE
    dCb_dt = -1 / V * mp / rho_p * 3 / R * De * dCp_drR;

    % Combine into a single vector
    dydt = [dCp_dt; dCb_dt];
end

[Torsten]

If kR*Cp in your document must be replaced by kR*Cp.^n for an n-th order reaction, the changes you made should be correct.

For the discretization, see Chapter 3 of the enclosed document.

[Aditya]

Thanks

Answer 2

[Torsten]

tau is a vector of values in your code:

tau = linspace(0, tau_end, 100); % Non-dimensional time

You can't use it as a symbolic variable afterwards without redefining it as symbolic:

syms tau Cbbar(tau)

But note that you now have overwritten the original tau.

Further, replace

dCbbar_dtau = -alpha * (Bi * (Cbbar - Cpbar_surface') - 3 * phi^2 * Cpbar_surface');

by

dCbbar_dtau = diff(Cbbar,tau) == -alpha * (Bi * (Cbbar - Cpbar_surface') - 3 * phi^2 * Cpbar_surface');

to define your differential equation.

But this will only work for Cpbar_surface being a single value, not a vector of values that depend on time.

Thus you won't succeed to integrate

dCbbar_dtau = -alpha * (Bi * (Cbbar - Cpbar_surface') - 3 * phi^2 * Cpbar_surface');

with a time-dependent function Cpbar_surface(t) (in this case a vector of values) symbolically.

Use pde1dm if you want to solve a PDE with a coupled ODE:

https://uk.mathworks.com/matlabcentral/fileexchange/97437-pde1dm

Or discretize your PDE in space, add the ODE and solve the complete system of ordinary differential equations using "ode15s". Look up "method-of-lines" for more details.