Solving a system of nonlinear, discretized partial derivatives (mass balances) using FSOLVE

Dear MATLAB community,

I have a few questions regarding the topic that I hope some of you could help me on. The problem to solve is as follows:

For my research, I want to numerically analyze a system of stationary, 2D mass balances (electrochemistry) which is made up of partial derivatives and boundary conditions. I have by now found no built-in MATLAB commands are suitable (the format of the equations are not suitable for the PDE-tool, and pdepe is for transient systems only). Now I intend to discretize all equations using the command gradient, and solve the resulting system of algebraic, nonlinear equations using FSOLVE. Overall, 12 unknown variables (arrays) are to be solved by a system of 21 equations & constraints. I have embedded the relevant pieces of code.

My questions are:

- Whenever I want to run the m.file for testing, the following error shows up: Undefined function or variable 'x'. Why? I've seen numerous examples on MATLAB where multiple variables are defined as a = x(1), b = x(2) after which they are used in an equation F. This is the code:

function EQNS = nernst_planck_ae(x)
% Define independent variables in terms of x
c = x(1);
d = x(2);
cTm = x(3);
etc. until x(12).

- In essence, I want the solver to give arrays as an output (the variable concentration along the x-axis amongst others). However, when giving the equations, which are also array valued, the following error introduces itself: In an assignment A(:) = B, the number of elements in A and B must be the same. I understand this, but how can I work around this? The relevant code of line is:

EQNS(4) = - J_ch - D_m * cTm(end) .* gradient(phi_iem)

- Although the command gradient allows me to take the central difference of the derivative for any array/matrix, I also want to implement some boundary constraints for the solver. More specifically, I want for f(x), df/fx at x = 0 to be 0. How do I implement such a constraint in the equations? Is the following line of code in the correct format?:

% Neumman conditions for the dcdx & J_ions
  dcdx = gradient(c);
  dddx = gradient(d);

    EQNS(10) = dcdx(1) == 0;                             
    EQNS(11) = dcdx(end) == -J_ions/(2*D_s);            
    EQNS(12) = dddx(1) == 0;                           
    EQNS(13) = dddx(end) == -J_ions/(2*D_s);

The EQNS are objective functions for fsolve.

I realize this is quite some text, though I tried to ask as concise as possible without letting out relevant information. I can upload the entire script if that helps. Also, I can reference the academic paper if needed.

In any case, thanks a lot for the effort in advance,

Jair

Answer 1

[Torsten]

First question:

From what you post we can not tell the reason for the error message. Maybe your call to "fsolve" will tell us more.

Second question:

EQN(4) is a scalar whereas - J_ch - D_m * cTm(end) .* gradient(phi_iem) seems to be a vector. So setting them equal is not possible.

You will have to reserve one equation for each component of the vector.

Third question:

    EQNS(10) = dcdx(1);                            
    EQNS(11) = dcdx(end) - (-J_ions/(2*D_s));            
    EQNS(12) = dddx(1) ;                           
    EQNS(13) = dddx(end) - (-J_ions/(2*D_s));

Best wishes

Torsten.

[Jair Dan]

Thanks for the quick answer!

1. The error display from the solver is as follows:

Index exceeds matrix dimensions.

Error in nernst_planck_ae (line 9)
cTm = x(3);

Error in fsolve (line 230)
            fuser =
            feval(funfcn{3},x,varargin{:});

Error in sol_nernst_planck (line 5)
solution = fsolve(EQNS,x0);

Caused by:
    Failure in initial objective function
    evaluation. FSOLVE cannot continue.

What suprises me is that the error only seems to come up for x(3) and beyond, and not for x(1) and x(2)?

2. I will first compute the vectors for the gradients, and then run a for loop to reserve an equation for each element of the gradient vector. Considering I have to that with the majority of the functions, there is no limit to the equations for fsolve?

3. Thanks, makes sense indeed.

[Torsten]

1.

What is the length of x0 ? If it's 2, then you know the reason why MATLAB errors.

2.

No limit, if you have enough time and RAM.

Best wishes

Torsten.

[Jair Dan]

I found out now length(x0) = 58, which is because I have inserted arrays as elements for the x0 list (so arrays as initial guesses next to scalars for some functions). I reckon that's not the way to go after reading to above.

I tried applying arrayfun(fun,x) for the EQNS, but the output would still be an array. Seems I have to run 'FSOLVE' through 'array fun' to get an array of results?

If there is some extra documentation on FSOLVE beyond the MATLAB documentation that you know of perhaps? I would be like to read further in it to understand the options with the code. Right now the challenge seems to use fsolve, while having arrays as initial guesses and outputs. The problem is also that not all arrays are of the same length, so I cannot use them as input for arrayfun or a loop simultaneously, yet the system should be solved as a whole together.

Thanks for your time and answers Torsten!

[Jair Dan]

This is what I had coded now for example:

EQNS(6) = - cTm(1) + sqrt(X^2+(2*c(end))^2);                % eq. 15 (c)
EQNS(7) = - cTm(end) + sqrt(X^2+(2*d(end))^2);              % eq. 15 (d)
EQNS(8) = ...
    arrayfun(@(phi_iem,cTm) gradient(phi_iem)...
    .*gradient(cTm)+ gradient(gradient(phi_iem)) .* ...
    cTm, phi_iem, cTm);                                     % eq. 11= eq.17
EQNS(9) = ...
    arrayfun(@(cTm,phi_iem) gradient(gradient(cTm))...
    - X.* gradient(gradient(phi_iem)), cTm, phi_iem);       % eq. 12=eq. 17

I know arrayfun outputs an array as well, so it's not solved yet. But I cannot run solve through arrayfun or a loop since the length is not the same for all arrays (length(c) = 10, length(phi_iem) = 6). How would I go about solving the system simultaneously then?

Also, somehow the solver still doesn't accept the x-variables to be solved:

function EQNS = ED_model(x)
%%Model parameters
% Define independent variables in terms of x
c = x(1);
d = x(2);
cTm = x(3);
phi_iem = x(4);
phi_d = x(5);
phi_c = x(6);
don_c = x(7);
don_d = x(8);
c_eff = x(9);
d_eff = x(10);
J_ions = x(11);
J_ch = x(12);

Error in fsolve (line 230)
            fuser =
            feval(funfcn{3},x,varargin{:});

Thanks again,

Jair

[Torsten]

count1 = length(array1)
for i=1:count1
  EQNS(8+i)=array1(i);
end

count2 = length(array2)
for i=1:count2
  EQNS(8+count1+i)=array2(i)
end 

...

Best wishes

Torsten.

[Jair Dan]

Thanks that helped a lot, I've got the system running (at least). Now I just have to keep finetuning it.