how Can I have a stable and converged solution using fsolve?
9 次查看(过去 30 天)
显示 更早的评论
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]);
0 个评论
回答(1 个)
Roger Stafford
2015-4-16
编辑:Roger Stafford
2015-4-16
"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.
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!