how Can I have a stable and converged solution using fsolve?

Hello everyone,

I have the following parameters and functions,I used fsolve to iteratively solve the nonlinear function. there is only one variable iteratively solved. my problem is that when I change the initial guess my solution (Hw) changes and I do not want it because solution must not depend on the initial guess. lets say, when I have initial guess of 0.9 or 0.8, my result is very close to the reality while when it is less than 0.5, it is unphisically meaningless. Is there any other solvers that I can use to have unique solution?

D=0.030; A=pi*(D^2)/4; ro_o=910; ro_w=1000; mu_w=0.001; Jo=[0.818929836 0.818929836 0.818929836 0.818929836 1.094285163 1.094285163 1.094285163 1.094285163 1.36964049 1.36964049 1.36964049 1.36964049 1.644995817 1.644995817 1.644995817 1.644995817];

Jw=[1.192024677 1.582473576 1.974263581 2.362342154 1.198574263 1.590988038 1.974419523 2.359348057 1.184083303 1.576578948 1.97928493 2.340822085 1.191837546 1.581070093 1.97398572 2.359226422];

mu_o=[0.677160616 0.677160616 0.677160616 0.677160616 0.677160616 0.677160616 0.677160616 0.677160616 0.677160616 0.677160616 0.677160616 0.677160616 0.677160616 0.677160616 0.677160616 0.677160616];

dpdx_exp=[1.246991 1.714304 2.17153 3.132296 1.764007 2.592897 2.985523 3.272178 2.384805 2.617821 2.993001 3.409082 2.267615 2.869545 3.286467 3.909495];

for k=1:16

Rew(k)=(((Jw(k)*D*ro_w)/mu_w)^-.2); sw=pi*D; B(k)=0.023*ro_w*Rew(k)*sw;

options = optimoptions('fsolve','Display','iter');

Hw(k) = fsolve(@(Hw) (((0.5*ro_o*(-3.6*log(6.9/(((D-(2*sqrt(Hw*A/pi)))*(Jo(k)/(1-Hw))*ro_o)/mu_o(k))))^-2.0)*(((Jo(k)/(1-Hw))-(Jw(k)/Hw))*((Jo(k)/(1-Hw))-(Jw(k)/Hw))))*(pi*(D-(2*sqrt(Hw*A/pi)))*(1/(1-Hw))))-((B(k)*(Jw(k)/Hw)^2)), 0.9, options)

Do(k)=D-(2*sqrt(Hw(k)*A/pi)); nu_o(k)=mu_o(k)/ro_o; Reo(k)=(Do(k)*(Jo(k)/(1-Hw(k))))/nu_o(k);

fow(k)=(-3.6*log(6.9/Reo(k)))^-2; tau_i(k)=0.5*ro_o*fow(k)*((Jo(k)/(1-Hw(k)))-(Jw(k)/Hw(k)))*(((Jo(k)/(1-Hw(k)))-(Jw(k)/Hw(k)))); si(k)=pi*Do(k); dpdx(k)=(tau_i(k)*si(k)/(A*(1-Hw(k))))/1000;

end

plot(dpdx,dpdx_exp,'o') hold plot([0 5],[0 5]);

Answer 1

[Roger Stafford]

"when I change the initial guess my solution (Hw) changes" This is a feature of 'fsolve' you will have to accept, Parham. Its solution will often depend on the initial estimate when there are multiple solutions. The remedy is to either use an analytic solution as obtained from 'solve', or else provide multiple initial estimates to 'fsolve' from which you can select the appropriate one.

Consider the very simple problem x^2-x-6 = 0. It has two solutions. If you give an initial estimate greater than .5, it will undoubtedly converge to the solution x = 3, but if the initial estimate is less than .5, it will converge to the solution x = -2. There is nothing you can do about this if you use iterative routines like 'fzero' or 'fxolve' except to use multiple initial estimates or else use some non-iterative method which provides multiple solutions.

One technique that is useful in cases such as yours with only one independent variable is to make a plot of the expression which you would like to have equal to zero, as a function of the independent variable. From the plot you can usually devise initial estimates that will behave as you wish in the presence of multiple zero crossings for providing accurate solutions.

[Parham Babakhani Dehkordi]

Thank you for you response,Roger. it helps a lot. I have another question. Is it possible to separate all the terms that contain Hw to prevent the complication of the above function? I mean, is there any way not to expilicitly put all therms in one equation and then pass the fsolve to it?