Why is ODE15s ignoring my boundary conditions?

Hi all, I'm attempting to solve the Fishers PDE by spatially discretising and then feeding the time derivative to ODE15s (therefore creating a system of ODEs). Here's the equation in full;

.

I've set the initial conditions as and then I have set the second element of this to so it respects the boundary conditions at . Here's my code:

space_steps=500;
T=100;
x = linspace(0,1,space_steps).';
T_span=[0;T]; %T_vector
n_0 = 0.6*(1-x.^2);
n_0(2)=0.6;

[T_out,Y_out] = ode15s(@(t,y) Fisher(t,y,space_steps),T_span,n_0);

Here's the function 'Fisher' that does all the work:

function output = Fisher(~,input,space_steps)

n = [input(2);input(2:end-1);0];    %dndx = 0 at x=0 so input(1)=input(2), n=0 at the x=1.
dx=1/(space_steps-1);               %Delta x

np = [n(2:end);NaN]; %n_{i+1}
nm=[NaN;n(1:end-1)]; %n_{i-1}

%Second order derivative with a forward and backward finite difference quotient at each boundary respectively
d2ndx2 = [(2*n(1)-5*n(2)+4*n(3)-n(4))/dx; np(2:end-1)-2*n(2:end-1)+nm(2:end-1);(2*n(end)-5*n(end-1)+4*n(end-2)-n(end-3))/dx]/dx^2;


dndt =  d2ndx2 + n.*(1-n); %sets up the equation 
output = dndt; %gives it to ODE15s
end 

Now for some reason the function or ODE15s is completely ignoring the two boundary conditions I am prescribing it at each run at the begining. Could anyone explain why this is and provide a suitable remedy? Thanks

[darova]

It's because of these guys

Ordinary means that equations have one independet variable (time maybe or something else) LINK

You have two independent! That means your equations is not ordinary. Your equations is partial

You can use pdetoolbox or try to solve it using FDM (Finite Difference Scheme)

[John D'Errico]

@darova - Actually, the spatial discretization will work, if done properly. The idea is to treat it as SYSTEM of coupled ODEs, where the inkage comes from the discretization in x. Thus you replace the partial derivatives in x with a finite difference. Then just asssume you have n_x variables to solve for at each time step, so a system of n_x ODES as a function of time.

I''ll try to offer an answer in a bit, if nobody else has done so already. That will mean verifying the discretization has been done correctly, and the boundary values treated properly (which seems to be the issue at hand.)

[Jacob Jepson]

@darova - thanks for your comment. However, if you read what I put, the spatial part is discretised which means the PDE actually gets transformed in to a system of ODEs, which ODE15s can solve. This is called the method of lines. As there becomes only one variable.

John, thank you but since I posted this question I have come up with a solution by taking advantage of the ‘conv’ function which produces the spatial discretisation by matrix form. Many thanks.

[Jacob Jepson]

@darova also might suggest reducing the amount italics and bolds in your comments in future so as not to come across as slightly condescending.

[darova]

Sorry, didn't notice that

[K Benis]

Hello,

I also faced the same problem in solving the follwing equations. I am wondering if you can explain me how can I solve this issue.

Independent variables are t and r. Dependent variables are CA (function of t), CAr (function of r,t).

There is a relation between q and CAr (q=f(CAr)). This relation can be a simple linear eqution (q=aCAr+b), or can be a complex equation like Eq. 8 or 9. Using this eqution the term q can be converted to CAr.

Using the method of lines, the derivatives in the PDE equation (Eq.3 ) are discretized into N+1 points on the spatial derivatives (r), where N is the number of grid points. I also used a second-order forward difference for B.C at r=0 (Eq. 5), and second-order backward difference at r=R (Eq. 6).

Below are the codes based on Eq. 7 (the simplest condition).

clear all
close all
clc
%
global DP DS KF RP RoP S V m EP nr r h

% other independent parameters
DP = 1.1e-10; % m2/s
DS = 5.17e-12; % m2/s
KF = 4.05e-5; % m/s
CA0 = 2e3; % g/m3
RP = (580/2)*10^-6; % m Particle radious
RoP = 957e3; % gr/m3
S = 6/(2*RP*RoP); % m2/g external surface area
V = 40e-5; % m3 L Volume of the solution
m = 20; % g added sorbent
EP = 0.352; % Porosity


% numerical calculation
nr = 15; % number of points in r direction
h = RP/(nr); % step size in r direction
t0 = 0; tf = 1200; %seconds
timerange =[t0 tf];
CAr0=zeros(nr+2,1);
CAr0(nr+2,1) = CA0;

for i=1:nr+1
    r(i) = h*(i-1);
end

%Linear
[t,CAr] = ode15s(@(t,CAr) F_linear(t,CAr), timerange, CAr0);

Function

function CArt = F_linear(t,CAr)

global KL DP DS KF RP RoP S V m EP nr r h

% B.C. at r=0 (center of particle) 2nd order forward difference
    CAr(1) =(4*CAr(2)-CAr(3))/3;
        
    for i=2:nr
    x1 = (CAr(i+1)-2*CAr(i)+CAr(i-1))/(h^2); %d2C/dr2
    x2 = (CAr(i+1)-CAr(i-1))/(2*h); %dc/dr
    x3 = (DP*(x1 + (2/(r(i))*x2))) + (DS*RoP*((0.124*x1)+ ((2*0.124)/(r(i))*x2)));
    x4 = EP + (0.124*RoP);
   
    CArt(i) = x3/x4; % Eq. 3
    end

   % B.C. at r=R (surface of particle) Calc CAr(n+1,t) 2nd order backward
    % ((3*CArnp1 - 4*CArn + 4*CArnm1)/(2*h)); %dc/dr
    slope = 0.124; % (Eq. 7)
    CArn = CAr(nr); CArnm1 = CAr(nr-1); CA = CAr(nr+2);
    y1 = ((2*DP)/h) + ((2*slope*RoP*DS)/h);
    y2 = (DP/(2*h)) + ((slope*RoP*DS)/(2*h));
    y3 = ((3*DP)/(2*h)) + ((3*slope*RoP*DS)/(2*h)) + KF;
    
    CArnp1 = ((KF*CA) + (CArn*y1) - (CArnm1*y2))/y3;
    
    CAr(nr+1) = CArnp1;
        % External mass transfer (Eq. 1)
    CArt(nr+2) = ((-m*KF*S)/V)*(CAr(nr+2)-CAr(nr+1));
    CArt = CArt(:);

end

Answer 1

[darova]

I rewrote your equations as following

I used this difference scheme

And here is the success i reached