MATLAB: Gradient algorithm in FSOLVE

boundary constraintderivativefsolvegradientpde

Hi all!
I have a question on the usage of the command/algorithm of gradient for equations in FSOLVE. For the systems of equations to be solved, I want to add boundary constraints on the derivatives of concentration c at x = 0 and x = end, so that dcdx(x=0) = 0 for example. For this purpose, I want to use the command gradient to compute discretized algebraic form of the derivative for numerical analysis. Concrete, my question is the following:
When pre-assigning the gradient of say 'c' to array 'dcdx', and then consequently using the array 'dcdx' in the equations for the objective function of FSOLVE, will it 'dcdx' pass on the algorithm of gradient or only the scalar values it calculated. Like in the following line of code (taken from the doc gradient):
dcdx(:,j) = 0.5*(c(:,j+1) - c(:,j-1));
The corresponding (piece of) code for my system is:
%%Model equations
% Discretize all balances into AEs, and re-arrange them to F(x) = 0
% Compute gradients for variables
[c_dx, c_dy] = gradient(c); 
[d_dx, d_dy] = gradient(d);
c_dx2 = gradient(c_dx); d_dx2 = gradient(d_dx);
phic_dx = gradient(phi_c); phid_dx = gradient(phi_d);
phic_dx2 = gradient(phic_dx); phid_dx2 = gradient(phid_dx);
phiem_dx = gradient(phi_iem); phiem_dx2 = gradient(phiem_dx);
cTm_dx = gradient(cTm); cTm_dx2 = gradient(cTm_dx);
% *** Algebraic equations ***
% Cell pair voltage & potentials
eq_16 = ...
    V_CP./V_T./2 - (phi_c(:,end) + phi_d(:,end) + don_c - ...
    don_d + phi_iem(:,end));                                % eq. 16 
% Donnan potentials
eq_14c = don_c - asinh(wX ./ (2.*c(:,end)));                % eq. 14
eq_14d = don_d - asinh(wX ./ (2.*d(:,end)));                % eq. 14 
    % Total ionic flux
  for j = 1:y_steps
      for i = 2:m_steps-1
          eq_11(j,i) = - J_ch(j) + (-D_m .* cTm(j,i) .* ...
              phiem_dx(j,i));                                % eq. 11
          eq_12(j,i) = - J_ions(j) + (-D_m .* ...
              (cTm_dx(j,i) - wX .* phiem_dx(j,i)));           % eq. 12
          eq_11_17(j,i) = phiem_dx(j,i) .* cTm_dx(j,i)...
              - (- cTm(j,i) .* phiem_dx2(j,i));               % eq. 11= eq.17
          eq_12_17(j,i) = cTm_dx2(j,i) - ...
              (wX .* phiem_dx2(j,i));                         % eq. 12= eq.17
      end
  end
% Membrane concentrations
eq_15c = - cTm(:,1) + sqrt(wX^2+(2*c(:,end)).^2);           % eq. 15 (c)
eq_15d = - cTm(:,end) + sqrt(wX^2+(2*d(:,end)).^2);         % eq. 15 (d)
% Neumman conditions for the dcdx & J_ions
bc_1 = c_dx(:,1);                                   % NC:dcdx(0,:) = 0                       
bc_2 = c_dx(:,end) - (-J_ions/(2*D_s));             % NC:dcdx(np,:)= eq.13  
bc_3 = d_dx(:,1);                                   % NC:dcdx(0,:)= 0
bc_4 = d_dx(:,end) + (-J_ions/(2*D_s));             % NC:dcdx(np,:)= eq.13
As you can see in the code, I first assign the gradients to arrays, and then use those arrays in the equations of the system. Will this work to find the right derivatives for the system using FSOLVE? Or should I use another way? If not, what other way is there?
Thanks a lot in advance,
Jair

Best Answer

Hello,
You have the right idea in my opinion. If you refer the documentation for fsolve
https://www.mathworks.com/help/optim/ug/fsolve.html#inputarg_fun
The function handle "fun", is a function that accepts a vector x and returns a vector F, the nonlinear equations evaluated at x.
So basically if you want to impose Boundary Conditions (BCs) that the solution x must satisfy, they will have to be returned as part of F.
gradient as such will only compute the numerical gradient at each grid point based on finite differences. It does not return a formula/function handle.
By specifying your BCs to be part of F,
bc_1 = c_dx(:,1);                                   % NC:dcdx(0,:) = 0                       
bc_2 = c_dx(:,end) - (-J_ions/(2*D_s));             % NC:dcdx(np,:)= eq.13  
bc_3 = d_dx(:,1);                                   % NC:dcdx(0,:)= 0
bc_4 = d_dx(:,end) + (-J_ions/(2*D_s));             % NC:dcdx(np,:)= eq.13     
F = [someConditions; bc_1 ; bc_2; bc_3; bc_4]
The solver will iterate through different values of x such that F(x) is zero or close to zero and when that happens, your BCs are automatically satisfied, and when that happens, the numerical values returned by gradient are the right ones.
Related Question