Gradient algorithm in FSOLVE
9 次查看(过去 30 天)
显示 更早的评论
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 个评论
采纳的回答
Mukul Rao
2017-4-24
Hello,
You have the right idea in my opinion. If you refer the documentation for fsolve
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.
更多回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Partial Differential Equation Toolbox 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!