Why is ODE15s ignoring my boundary conditions?

Question

Jacob Jepson 2020-4-16

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/518450-why-is-ode15s-ignoring-my-boundary-conditions

评论： K Benis 2020-6-9

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

6 个评论
显示 4更早的评论隐藏 4更早的评论

darova 2020-4-17

Sorry, didn't notice that

K Benis 2020-6-9

在 MATLAB Online 中打开

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