Solving system of coupled nonlinear discretized equations with fsolve using boundary conditions

Hi all,

I have an equation system of 5 coupled, nonlinear equations which are (spatially) discretized by first order central finite difference method (FDM) and I want to solve it with Matlab's fsolve function. Additionally, I need to implement boundary conditions.

The script for defining the fsolve function and the boundary condition is as follows:

% Parameterization of variables
Curr_dens=9000;
Fa=96485.33;
Z_bc_H2O=0.8;           %Initial/boundary molar concentration species 1
Z_bc_H2=1-Z_bc_H2O;     %Initial/boundary molar concentration species 2
Nz=20;              %number of grid points in the discretized mesh (Nz,1) mesh
dz=2.6316*1E-5;     %delta z of the grid
numb_eq=5;          %total number of equations to be solved
D_k_i=1E-5*[7.49;2.24];
D_B=1.3806*1E-4;
vars=[8.3145; 1073.15; 101325];


% Definition of the boundary conditions
N_bc_ij=zeros(Nz, 2);
N_bc_ij(end,1)=Curr_dens/(2*Fa);             %boundary condition for flux of species 1
N_bc_ij(end,2)=-Curr_dens/(2*Fa);            %boundary condition for flux of species 2
Z_bc_ij=zeros(Nz, 2);
Z_bc_ij(1,1)=Z_bc_H2O;          %Boundary condition for molar concentration of species 1 
Z_bc_ij(1,2)=Z_bc_H2;           %Boundary condition for molar concentration of species 2

% Definition of the solution function
fun=@(f) mymodel (f, Nz, dz, numb_eq, vars, ...
        D_k_i, D_B, N_bc_ij, Z_bc_ij);

%Providing initial guesses
f_guess=zeros(Nz,numb_eq);
f_guess(end-1:end,1)=Curr_dens/(2*Fa);       %Molar flux species 1    
f_guess(end-1:end,2)=-Curr_dens/(2*Fa);      %Molar flux species 2
f_guess(1:2,3)=Z_bc_H2O;                    %molar concentration species 1                
f_guess(1:2,4)=Z_bc_H2;                     %molar concentration species 2
f_guess(:,5)=1;                             %sum of molar concentrations must be 1
f_guess=f_guess(:);

% Solution of the equation system with fsolve
options=optimoptions('fsolve', 'Display','iter', 'FunctionTolerance',1e-5, ...
    'MaxIterations', 1000, 'MaxFunctionEvaluations', 1e7);
f_solution=fsolve(fun,f_guess, options);

%Reshaping of the solution
N_H2O=f_solution(1 : 1*Nz).';
N_H2=f_solution(1*Nz+1 : 2*Nz).';
Z_H2O=f_solution(2*Nz+1 : 3*Nz).';
Z_H2=f_solution(3*Nz+1 : 4*Nz).';
Z_total=f_solution(4*Nz+1 : end).';

and the fsolve function block with my equations (to calculate molar flux and molar fraction of two species) as follows:

function [F] = mymodel (f, Nz, dz, numb_eq, vars, ...
        D_k_i, D_B, N_bc_ij, Z_bc_ij)

% Reshaping of the solution vector into the vectors for each of the species
N_H2O=f(1 : 1*Nz);                %Solution vector for the flux of species 1           
N_H2=f(1*Nz+1 : 2*Nz);            %Solution vector for the flux of species 2
Z_H2O=f(2*Nz+1 : 3*Nz);           %Solution vector for the molar concentration of species 1
Z_H2=f(3*Nz+1 : 4*Nz);            %Solution vector for the molar concentration of species 2
Z_total=f(4*Nz+1 : 5*Nz);         %Solution vector for the total molar concentration
% Creation of the matrix containing the 5 equations over the number of grid poins
F=zeros(Nz,numb_eq);               

for z=2:Nz-1            
%from the second to the second last cell because of central FD discretization
    % Flux equations for the both species
    F(z,1)=(N_H2O(z+1)-N_H2O(z-1))/(2*dz);           
    F(z,2)=(N_H2(z+1)-N_H2(z-1))/(2*dz);

    %Implicit expressions for fluxes and molar concentrations for the both
    %species
    F(z,3)=(Z_H2O(z+1)-Z_H2O(z-1))/(2*dz) + vars(1)*vars(2)/vars(3) * (...
        (Z_H2(z)*N_H2O(z) - Z_H2O(z)*N_H2(z))/D_B + N_H2O(z)/D_k_i(1));
    F(z,4)=(Z_H2(z+1)-Z_H2(z-1))/(2*dz) + vars(1)*vars(2)/vars(3) * (...
        (Z_H2O(z)*N_H2(z) - Z_H2(z)*N_H2O(z))/D_B + N_H2(z)/D_k_i(2));

    %For assuring that the total concentration is 1
    F(z,5)=Z_total(z)-(Z_H2O(z)+Z_H2(z));
end

%Enforcing boundary conditions for the fluxes for both species  
F(end,1)=N_H2O(end)-N_bc_ij(end,1);
F(end,2)=N_H2(end)-N_bc_ij(end,2);
F(1,3)=Z_H2O(1)-Z_bc_ij(1,1);  
F(1,4)=Z_H2(1)-Z_bc_ij(1,2);
F(:,5)=Z_total(:)-ones(Nz,1);

F=F(:); %Reshaping of the solution vector for fsolve
end

The main problem I face is that fsolve sticks to the initial guess and almost completely ignores any boundary conditions I try to implement. Changes in the calculation methodology or in the equations also only lead to a very small deviation of the solution from the initial guess.

Can anyone please help me with this problem?Is the implementation with the fsolve function of Matlab a good way or would another Matlab solver be more suitable?

Thank you very much in advance.

Answer 1

[Torsten]

Use "bvp4c".

Additionally, I don't understand why you set

F(z,1)=(N_H2O(z+1)-N_H2O(z-1))/(2*dz);           
F(z,2)=(N_H2(z+1)-N_H2(z-1))/(2*dz);

This is a discretization of

dN_H2O/dz = 0
dN_H2/dz = 0

giving

N_H2O(z) = N_H2O(z=0)
N_H2(z) = N_H2(z=0)

Is that really what you want ?

[Stefan]

Dear Torsten,thank you very much for your reply. The problem I am facing with bvp4c is, that I only have boundary conditions at the interval end for the molar flux (equations 1 and 2) and boundary conditions at the begin of the interval for the molar concentration (equations 3 and 4). So I do not have boundary conditions for the given equations on both of the interval. Is the bvp4c solver still applicable under these circumstances, since according to my information this solver requires boundary conditions for each equation on both sided of the intervall?

Regarding the discretization of the flux: yes, for now this is what I need, there is one source term at the end of the interval and the fluxes do not change over the interval

[Torsten]

If you include your equations with boundary conditions in a mathematical notation, I can tell you whether bvp4c is applicable in your case.

If you have first-order odes (as seems to be the case), it's mathematically impossible to prescribe the flux. This only works for second-order odes.

[Stefan]

These are the boundary conditions I have:

(index1=H2O and index2=H2 in the previous chat)

For better understanding, these are the equations (system) I want to solve for i=1...2:

According to your suggestion of mathematical prescription of the molar flux for second-order odes: I could try to rearrange the second equation for Ni and substitute it into the first equation, or would another approach be more effective in your opinion?

[Torsten]

N1(1) = 0.05;
N1(2) = -0.05;
P = 101325;
R = 8.3145;
T = 1073.15;
D_k_i(1) = 7.49e-5;
D_k_i(2) = 2.24e-5;
D_B = 1.3806e-4;
fun_ode45 = @(z,Z) [-R*T/P*((Z(2)*N1(1)-Z(1)*N1(2))/D_B + N1(1)/D_k_i(1));...
                    -R*T/P*((Z(1)*N1(2)-Z(2)*N1(1))/D_B + N1(2)/D_k_i(2))];
Nz=20;              
dz=2.6316*1E-5;     
zspan = [ 0 Nz*dz ];
Z0 = [0.8;0.2];
[z,Z] = ode45(fun_ode45,zspan,Z0);
figure(1)
plot(z,Z)
grid on

zmesh = 0:dz:Nz*dz;
profun = @(z,N,Z)[0;0;...
                  -R*T/P*((Z(2)*N(1)-Z(1)*N(2))/D_B + N(1)/D_k_i(1));...
                  -R*T/P*((Z(1)*N(2)-Z(2)*N(1))/D_B + N(2)/D_k_i(2))];
fun_bvp4c = @(z,u)profun(z,u(1:2),u(3:4));
bcfun = @(Nl,Zl,Nr,Zr)[Nr(1)-N1(1);Nr(2)-N1(2);Zl(1)-Z0(1);Zl(2)-Z0(2)];
bc_bvp4c = @(ul,ur)bcfun(ul(1:2),ul(3:4),ur(1:2),ur(3:4));
solinit = bvpinit(zmesh, [N1(1); N1(2); Z0(1); Z0(2)]);
sol = bvp4c(fun_bvp4c,bc_bvp4c,solinit);
figure(2)
plot(sol.x,[sol.y(3,:);sol.y(4,:)])
grid on

[Stefan]

Thank you so much for your answers, you really helped me a lot.