Initial conditions for PDEPE

Hi community,

I am trying to solve 9 PDEs using pdepe routine for which the initial condition for all these equations is zero, Ci(t=0) = 0. My solution mesh is as:

% Solution mesh
x = linspace(0,1,50); 		% space
t = linspace(0,100,100);	% time

I attempted to code the initial conditions as:

% icfun_sr() defines the initial conditions
function u0 = icfun_sr(x)
u0 = zeros(9,1);
end

which I can then pass to solve the equations as:

% Solve equations
m=0;
sol = pdepe(m, @pdefun_sr, @icfun_sr, @bcfun_sr, x, t, [], A, B, D);

The error thrown as

Error using sr_code>icfun_sr

Too many input arguments.

Error in pdepe (line 229)

temp = feval(ic,xmesh(1),varargin{:});

Error in sr_code (line 32)

sol = pdepe(m, @pdefun_sr, @icfun_sr, @bcfun_sr, x, t, [], A, B, D);

Full code below:

close all
clear all
clc

global A B D Ncomp

% Constants
Rtime   =  2;                  
u     =  0.7;                             
L     =  1;                 
D_ax  = 1.5e-3;              
Bo    = u*L/D_ax;            
density = 2000;                
poros = 0.45;
Ncomp = 9;

% mesh
x = linspace(0,1,50); 		% space
t = linspace(0,100,100);	% time

% Other parameters;

A = -1/Rtime;
B = 1/Rtime*Bo;
D = density/poros;

% Solve equation
m=0;
pdef = @(x, t, u, dudx) pdefun_sr(x, t, u, dudx);
sol = pdepe(m, pdef, @icfun_sr, @bcfun_sr, x, t);

% Extract solutions
for i = 1:Ncomp
	u(i) = sol(:,:,i);
end

waterfall(x,t,u(1)), map=[0 0 0]; colormap(map), view(15,60)


function [c,f,s] = pdefun_sr(x, t, u, dudx)

global A B D Ncomp

c = ones(Ncomp,1); 	% diagonal terms

k = ones(12,1);
u(10) = 6; % or u(10) = 6 - exp(-10*t);

f = 	[0;		% diffusion term                
 	    0;                                      
	    0;                                      
	    0;                                      
	    A*u(5) + B*dudx(5);                     
	    0;                                      
	    A*u(7) + B*dudx(7);                     
	    A*u(8) + B*dudx(8);                     
	    A*u(9) + B*dudx(9)                     
	    ];

% source term
s = 	[k(7).*u(3) - k(1).*u(1); 
	    k(9).*u(6)   -  k(2).*u(2); 
	    k(12).*u(10) - (k(3) + k(5) + k(8) + k(7)).*u(3); % fails at u(10)
	    k(10).*u(6)  +  k(11).*u(5)  - k(6).*u(4) - k(4).*u(4); 
	    D*( k(6).*u(4) - k(11).*u(5) );
	    k(8).*u(3)   -  (k(9) + k(10)).*u(6);
	    D*( k(1).*u(1) );  
	    D*( k(2).*u(2) + k(3).*u(3) );
	    D*( k(5).*u(3) + k(4).*u(4) )
	    ];
end

% icfun_sr() defines the initial conditions
function u0 = icfun_sr(x)

global Ncomp

u0 = zeros(Ncomp,1);
end

function [pl,ql,pr,qr]  = bcfun_sr(xl,ul,xr,ur,t,A,B,D)
% bcfun_sr() defines the initial conditions

global Ncomp

pl = zeros(Ncomp,1).*ul;
ql = zeros(Ncomp,1);
pr = zeros(Ncomp,1);
qr = ones(Ncomp,1); % or ones(dudx)*Ncomp times
end

Answer 1

[Torsten]

Maybe you could add a little diffusivity to the components with spurious oscillations (those with f=0). If this is not possible: I gave you two alternatives.

global A B D Ncomp
% Constants
Rtime   =  2;                  
u     =  0.7;                             
L     =  1;                 
D_ax  = 1.5e-3;              
Bo    = u*L/D_ax;            
density = 2000;                
poros = 0.45;
Ncomp = 9;
% mesh
x = linspace(0,1,500); 		% space
t = linspace(0,100,100);	% time
% Other parameters;
A = -1/Rtime;
B = 1/Rtime*Bo;
D = density/poros;
% Solve equation
m=0;
sol = pdepe(m, @pdefun_sr, @icfun_sr, @bcfun_sr, x, t);
u = sol(40,:,6);
plot(x,u)

function [c,f,s] = pdefun_sr(x, t, u, dudx)
global A B D Ncomp
c = ones(Ncomp,1); 	% diagonal terms
k = ones(12,1);
u10 = 6; % or u10 = 6 - exp(-10*t);
f = 	[0;		% diffusion term                
    0;                                      
    0;                                      
    0;                                      
    B*dudx(5);                     
    0;                                      
    B*dudx(7);                     
    B*dudx(8);                     
    B*dudx(9)                     
    ];
% source term
s = 	[k(7).*u(3) - k(1).*u(1); 
    k(9).*u(6)   -  k(2).*u(2); 
    k(12).*u10 - (k(3) + k(5) + k(8) + k(7)).*u(3); % fails at u(10)
    k(10).*u(6)  +  k(11).*u(5)  - k(6).*u(4) - k(4).*u(4); 
    A*dudx(5) + D*( k(6).*u(4) - k(11).*u(5) );
    k(8).*u(3)   -  (k(9) + k(10)).*u(6);
    A*dudx(7) + D*( k(1).*u(1) );  
    A*dudx(8) + D*( k(2).*u(2) + k(3).*u(3) );
    A*dudx(9) + D*( k(5).*u(3) + k(4).*u(4) )
    ];
end
% icfun_sr() defines the initial conditions
function u0 = icfun_sr(x)
global Ncomp
u0 = zeros(Ncomp,1);
end
function [pl,ql,pr,qr]  = bcfun_sr(xl,ul,xr,ur,t,A,B,D)
% bcfun_sr() defines the initial conditions
global Ncomp
pl = ul;
ql = zeros(Ncomp,1);
pr = zeros(Ncomp,1);
qr = ones(Ncomp,1); % or ones(dudx)*Ncomp times
end

Answer 2

[Bill Greene]

The problem is that you are using an old, undocumented form of the of the pdepe function

to pass additional parameters to the functions called by pdepe. If you use this form, you must

include the extra parameters in ALL of the called functions, e.g.

function u0 = icfun_sr(x,A,B,C)

The reason that this old form of pdepe is undocumented is that anonymous functions provide a

superior way of passing extra parameters to called functions. In your example, you could define

the anonymous function

pdef = @(x, t, u, dudx) pdefun_sr(x, t, u, dudx, A, B, D);

and call pdepe

sol = pdepe(m, pdef, @icfun_sr, @bcfun_sr, x, t);

But even after changing your code to either of the two solutions, you will still get errors due

to the variable k being undefined in function pdefun_sr

.

[Matterazzi]

OK, thank you very much @Bill Greene for identifying this problem. I will continue working on the code and define all unknowns.

[Torsten]

You will get problems with the components of your solution vector for which you prescribed f = 0.

For the other components, take the convective parts of the fluxes into the source term s and don't leave them in f. The reason is that your formulation from above will collide with the boundary condition at x=xr.

[Matterazzi]

@Torsten, a problem manifests as you anticipated.

Error using pdepe (line 293)
Spatial discretization has failed. Discretization supports only parabolic and elliptic
equations, with flux term involving spatial derivative.

Error in sr_code3 (line 32)
sol = pdepe(m, pdef, @icfun_sr, @bcfun_sr, x, t);

According to your suggestion should the new f and s then be formulated in the following way?

f =    [0		% diffusion term
        0 
        0 
        0
        0
        0
        0
        0
        0
        ];


s = 	[k(7).*u(3) - k(1).*u(1)    % source term
        k(9).*u(6)   -  k(2).*u(2) 
        k(12).*u(10) - (k(3) + k(5) + k(8) + k(7)).*u(3)
        k(10).*u(6)  +  k(11).*u(5)  - k(6).*u(4) - k(4).*u(4) 
        A*u(5) + B*dudx(5) + D*( k(6).*u(4) - k(11).*u(5) )
        k(8).*u(3)   -  (k(9) + k(10)).*u(6)
        A*u(7) + B*dudx(7) + D*( k(1).*u(1) ); 
        A*u(8) + B*dudx(8) + D*( k(2).*u(2) + k(3).*u(3) )
        A*u(9) + B*dudx(9) + D*( k(5).*u(3) + k(4).*u(4) )
        ];
end

I am a beginner at this.

[Torsten]

You cannot use MATLAB's "pdepe" in its generic form.

Either you discretize your equations in space on your own or you try Bill Greene's code that can handle additional ODEs:

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

I think he will support you here if you have problems with his code.

And what I meant was to put the

A*dudx(5),A*dudx(7),A*dudx(8) and A*dudx(9) 

(the convective part) in s, not

B*d^2u/dx^2(5),B*d^2u/dx^2(7),B*d^2u/dx^2(8) and B*d^2u/dx^2(9) 

(the diffusive part).

[Bill Greene]

I suggest you post the equations you are trying to solve in mathematical form(you can use LaTex notation for this).

I'm fairly certain that your latest error is due to

pl = zeros(Ncomp,1).*ul;
ql = zeros(Ncomp,1);

not being a valid definition of boundary conditions at the left end.

[Matterazzi]

The nine equations are:

in the third equation can be taken from it's boundary conditions as

For all other components, the initial conditions and the boundary conditions and hold. Paramers A, B, and D remain as in the code above.

Thank you very much!

[Torsten]

A boundary condition for the second derivative (d^2C/dx^2 = 0 at x = 1) is impossible to prescribe. This follows from the theory of partial differential equations.

[Bill Greene]

I assumed that the error Torsten pointed out in your mathematical description was just a typo on your part-- that you really meant zero-flux at x=1. I also assumed that your left-end BC was also a typo-- that you really want C_i=0 at x=0. Accordingly I changed the pl definition in your BC function to

pl = ul;

This allows pdepe to complete the analysis. However C_1, C_2, C_3, C_4, and C_6 are showing spurious oscillations. This is likely due to the lack of a diffusion term in these equations.

[Matterazzi]

I have tried the @Bill Greene's pde1dm code with the problem and also used @Torsten's modification of the code and I get worthwhile results. This is what I exactly wanted. Many many thanks to you coders. You are really good at what you are doing! :)

[Torsten]

Your value for B is approximately 233 ! I expected something in the order of D_ax. Don't you think this must be wrong ?

Setting

f = 	[B*1e-7*dudx(1);		% diffusion term                
    B*1e-7*dudx(2);                                      
    B*1e-7*dudx(3);                                      
    B*1e-7*dudx(4);                                      
    B*dudx(5);                     
    B*1e-7*dudx(6);                                      
    B*dudx(7);                     
    B*dudx(8);                     
    B*dudx(9)                     
    ];

instead of

f = 	[0;		% diffusion term                
    0;                                      
    0;                                      
    0;                                      
    B*dudx(5);                     
    0;                                      
    B*dudx(7);                     
    B*dudx(8);                     
    B*dudx(9)                     
    ];

introduces a slight amount of diffusivity that prevents the spurious oscillations to happen (at least for the value of B you use at the moment).

[Matterazzi]

@Torsten to be honest, the values that go into B are some unfounded estimates that I came up with to focus on the maths. These will be fine tuned. It is a very brilliant idea of introducing infinitesimal diffusivity to deter spurious oscillations in those equations where f = 0. Thanks a bunch.

[Torsten]

Just out of interest: Why do components 5,7,8 and 9 move with velocity A and diffuse while the other components are spatially fixed and don't diffuse ?

[Matterazzi]

They are different phases. Components 1, 2, 3, 4, and 6 are fixed to certain positions and have lost one of the degrees of translational freedom. Components 5, 7, 8, and 9 interact within the system by virtue of their free nature, and retain their translational energy hence the convective, diffusive, and source terms all appear in their modeling.

[Matterazzi]

@Torsten, I am not sure if I should formulate this as a different question but I will ask here since it is a continuation of the problem above. Why is k equal to initial guess k0? I have tried several optimization approaches but nothing seems to solve the problem.

close all
clear
clc

% This problem does not solve because the function outpu does not change at all
% when your optimization parameters change a tiny amount. 

% Experimental data

Y=load('sr_data.txt'); % t C16O18O  H218O   CH3C16O18OH     C18O2
tdata    = Y(1:end,1);         
c9       = Y(1:end,2);         
c8       = Y(1:end,3);        
c5       = Y(1:end,4);        
c7       = Y(1:end,5);        
c        = [c5,c7,c8,c9];         % 36x4


k0 = [.02; .03; .01; .06; .03; .002; .4; .5; .01; .05; .1; .1];


% Three options to minimize
% lsqcurvefit returns error: Initial point is a local minimum
% lsqnonlin returns error: 'Exiting due to infeasibility:  11 lower bounds exceed the corresponding upper bounds.'
% patternsearch: from the web problem???

options = optimoptions('lsqcurvefit','FiniteDifferenceStepSize',5e-3, 'Algorithm','trust-region-reflective');
lb = zeros(12,1);
%lb = [];
ub = [];
[k,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,k0,tdata,c,lb,ub,options)

Initial point is a local minimum.

Optimization completed because the size of the gradient at the initial point 
is less than the value of the optimality tolerance.

k = 12×1

0.0200 0.0300 0.0100 0.0600 0.0300 0.0020 0.4000 0.5000 0.0100 0.0500

Rsdnrm = 9.5548e+04

Rsd = 36×4

-0.0024 -0.0030 -0.0046 0 -0.0058 0.0433 3.7644 3.7552 -0.0061 0.2084 7.9363 7.9050 -0.0051 0.4801 11.6359 11.5685 -0.0052 0.8454 14.9631 14.8392 -0.0179 1.2656 17.8793 17.7765 -0.0336 1.7379 20.5297 20.4439 -0.0593 2.2617 22.9251 22.7920 -0.1002 2.8260 25.0457 24.9382 -0.1277 3.4086 26.9769 26.8719

ExFlg = 1

OptmInfo = struct with fields:
      firstorderopt: 0
         iterations: 0
          funcCount: 13
       cgiterations: 0
          algorithm: 'trust-region-reflective'
           stepsize: 1
            message: 'Initial point is a local minimum.↵↵Optimization completed because the size of the gradient at the initial point ↵is less than the value of the optimality tolerance.↵↵<stopping criteria details>↵↵Optimization completed: The final point is the initial point.↵The first-order optimality measure, 0.000000e+00, is less than↵options.OptimalityTolerance = 1.000000e-06.'
       bestfeasible: []
    constrviolation: []

Lmda = struct with fields:
         lower: [12×1 double]
         upper: [12×1 double]
         eqlin: []
       ineqlin: []
      eqnonlin: []
    ineqnonlin: []

Jmat = 
   All zero sparse: 144×12

Cfit = kinetics(k, tdata);

% Plotting
plot(tdata, c, 'x')
hold on
plot(tdata, Cfit, '-',"LineWidth",1.2)

function C = kinetics(k,t)
global A B D Ncomp
% Constants
tau   =  11.8;                  
u       =  .9;                             
L       =  .3;                 
D_ax    = 1.5e-3;              
Bo      = 8.9;%u*L/D_ax;            
rho_B = 2000;                
eps_B   = 0.45;
Ncomp   = 9;

% solution mesh
x = linspace(0,1,20); 		% space
t = linspace(0,100,36);     % time

% Other parameters;
A = -1/tau;
B = 1/tau*Bo;
D = rho_B/eps_B;

% Solve equation
m=0;
options=odeset('RelTol',1e-4,'AbsTol',1e-4,'NormControl','off','InitialStep',1e-7);
sol = pdepe(m, @pdefun_sr, @icfun_sr, @bcfun_sr, x, t,options);

function [c,f,s] = pdefun_sr(x, t, u, dudx)
% global A B D Ncomp
c = ones(Ncomp,1);                  % diagonal terms
k = ones(12,1)*1E-2;
u10 = 6;
f = [B*1e-7*dudx(1);                % diffusion term            
    B*1e-7*dudx(2);                            
    B*1e-7*dudx(3);                            
    B*1e-7*dudx(4);                             
    B*dudx(5);           
    B*1e-7*dudx(6);                               
    B*dudx(7);           
    B*dudx(8);          
    B*dudx(9)            
    ];

s = 	[k(7).*u(3) - k(1).*u(1);   % source term
        k(9).*u(6)   -  k(2).*u(2); 
        k(12).*u10 - (k(3) + k(5) + k(8) + k(7)).*u(3); % fails at u(10)
        k(10).*u(6)  +  k(11).*u(5)  - k(6).*u(4) - k(4).*u(4); 
        A*dudx(5) + D*( k(6).*u(4) - k(11).*u(5) );
        k(8).*u(3)   -  (k(9) + k(10)).*u(6);
        A*dudx(7) + D*( k(1).*u(1) );  
        A*dudx(8) + D*( k(2).*u(2) + k(3).*u(3) );
        A*dudx(9) + D*( k(5).*u(3) + k(4).*u(4) )
        ];
end

% Initial conditions
function u0 = icfun_sr(x)

u0 = zeros(Ncomp,1);
end

% Boundary conditions
function [pl,ql,pr,qr]  = bcfun_sr(xl,ul,xr,ur,t,A,B,D)

pl = ul;
ql = zeros(Ncomp,1);
pr = zeros(Ncomp,1);
qr = ones(Ncomp,1); % or ones(dudx)*Ncomp times
end

% Extract modeled data used for fitting experimental data (at reactor exit)
u5 = sol(:,20,5);
u7 = sol(:,20,7);
u8 = sol(:,20,8);
u9 = sol(:,20,9);

C = [u5 u7 u8 u9]; % very good
end

[Torsten]

Why is k equal to initial guess k0?

Because you overwrite the parameters supplied by lsqcurvefit in this line:

k = ones(12,1)*1E-2;

[Matterazzi]

I am slapping myself now for being so blind. Thank you very much @Torsten.