Gradient algorithm in FSOLVE

Question

Jair Dan 2017-4-21

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/336473-gradient-algorithm-in-fsolve

评论： Jair Dan 2017-4-25

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Mukul Rao 2017-4-24

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/336473-gradient-algorithm-in-fsolve#answer_264306

在 MATLAB Online 中打开

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.

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Jair Dan 2017-4-25

Thanks for the answer Mukul!

I had followed the guidelines in your answer indeed, so I think the answer should be ok.

So if I understand correctly, my BC's are now rather values that are either 0 or not, which are iteratively changed by FSOLVE according to new values for x. In (simple) code:

bc_1 = scalar(x) or vector(x) == 0,

rather than bc_1 = formula(x) == 0?

请先登录，再进行评论。