How to set this time derivative boundary condition using MOL?

LHS = Left hand side aka

BC = Boundary condition

Hi,

I have a system of coupled raection diffusion PDEs as such.

For my species I have a LHS at such as

I am trying to solve this using the method of lines by using ODE15s solver for which I have written the function as such. But I do not know how to incorporate the term into the script and how to assign it to be my LHS boundary condition.

k1=0.073317274764052;
k2=1.052512778689219e-04;
tau=0.466500000000000;
tau_2=0.699700000000000; 
dchi=0.100000000000000;
N=11; 
tspan_2=linspace(tau, tau_2, N);
[t, CmBC] = ode15s(@(t, Cm2) BCodefun(t, Cm2, k1, k2, dchi), tspan_2, 0.9, []);
function dCdt2BC = BCodefun(t, CmBC, k1, k2, dchi)
         dCdt2BC  = (k1/k2).*(((CmBC).*(CmBC^2-CmBC))/dchi); 
end

This gives me a 1*11 matrix (code till above should function). You can see below how I am using a multi point approximation to set a no flux BC at my left hand. I am not sure how to assign my ODE output as my left hand side BC in my MOL implementation.

% What I have been using to set a no flux boundary condition on the LHS
Cm2(1)  = (-Cm2(3)+4*Cm2(2))./3;
% WHat I think I should do to set the output from my ODE solution to set my
% LHS BC
Cm2(1)  = CmBC;

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

[Torsten]

Use "pdepe" for your system if you are not experienced in programming numerical methods to solve PDEs.

[Hashim]

Hi Torsten,

Thank you for the reply. I have been using pdepe so far but as far as I can tell it is not possible or at least not the easiest to implement such boundary conditions which involve a time derivative in addition to the flux statement. It was actually one of your replies on one of the questions that prompted me to look into the method of lines. Here it is.

Solving PDE involving boundary condition with partial derivatives with respect to space and time. - MATLAB Answers - MATLAB Central (mathworks.com)

I would like to state here that I have replicated my previous simulations in MOL with ease so far.

Coule be that I have misunderstood your reply.

Thanks again

[Torsten]

Is c_m_ox the same as C_m_ox in your notation from above ?

If yes: What is the boundary condition for c_m_ox at 0 in its continuous form ? What you write above already seems to include a discretization in space.

[Hashim]

Yes, thank you for pointing it out.

It is as such

Yes sorry about that.

[Torsten]

Again a delta(Chi) appears in the denominator. This cannot be correct for a continuous form of the boundary condition.

[Hashim]

Yes corrected.

[Torsten]

Your solution variables within ode15s are C_s(1), C_s(2),...,C_s(N),C_m_ox(1), C_m_ox(2),...,C_m_ox(N) where N is the number of discretization points in space.

So what is the problem of implementing the ODE for C_m_ox in Chi=0 as

dC_m_ox/d(tau_M) = k1/k2 * (C_m_ox(2)-C_m_ox(1))/dChi * (C_m_ox(1)*(C_m_ox(1)-1))

?

[Hashim]

There is no in this boundary condition implementation only . Also, just to clarify I am using N to discretize both the spatial as well as the temporal domain.

This is what I am trying to do Cm2 being .

[t, Cm2(1)] = ode15s(@(t, Cm2) BCodefun(t, Cm2, k1, k2, dchi), tspan_2, 1, []);
function dCdt2BC = BCodefun(t, Cm2, k1, k2, dchi)
            dCdt2BC  = (k1/k2).*((Cm2-Cm2)/dchi)*(Cm2^2-Cm2); 
end

If I go your route I get the following error.

Index exceeds the number of array elements. Index must not exceed 1.

Error in PcWise_NDOCP_Pulse2v0/fpde2/BCodefun (line 133)
            dCdt2BC  = (k1/k2).*((Cm2(2)-Cm2(1))/dchi)*(Cm2(1)^2-Cm2(1));

Error in PcWise_NDOCP_Pulse2v0>@(t,Cm2)BCodefun(t,Cm2,k1,k2,dchi) (line 123)
    [t, Cm2(1)] = ode15s(@(t, Cm2) BCodefun(t, Cm2, k1, k2, dchi), tspan_2, 1, []);

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

Error in PcWise_NDOCP_Pulse2v0/fpde2 (line 123)
    [t, Cm2(1)] = ode15s(@(t, Cm2) BCodefun(t, Cm2, k1, k2, dchi), tspan_2, 1, []);

Error in PcWise_NDOCP_Pulse2v0>@(t,C2)fpde2(t,C2,alpha,kappa,eta,gamma,mu,k1,k2,dchi,N,tspan_2) (line 104)
[t, C2]      = ode15s(@(t, C2) fpde2(t, C2, alpha, kappa, eta, gamma, mu, k1, k2, dchi, N, tspan_2), tspan_2, IC2, []);

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

Error in PcWise_NDOCP_Pulse2v0 (line 104)
[t, C2]      = ode15s(@(t, C2) fpde2(t, C2, alpha, kappa, eta, gamma, mu, k1, k2, dchi, N, tspan_2), tspan_2, IC2, []);

Needless to say my implementation isn't working either.

[Torsten]

Reread how the method-of-lines works to solve partial differential equations.

As I wrote you don't solve for a single unknown value of the concentration field as you seem to try with the command

[t, Cm2(1)] = ode15s(@(t, Cm2) BCodefun(t, Cm2, k1, k2, dchi), tspan_2, 1, []);

but for the complete profile of C_m_ox and C_s over the spatial domain.

That's why I wrote you have 2*N solution variables where N is the number of spatial discretization points.

As far as I know, the following PDE solver similar to pdepe supports ODEs as boundary conditions:

https://uk.mathworks.com/matlabcentral/fileexchange/97437-pde1dm?tab=reviews

Maybe you want to test it for your problem.