Yet another solver that uses the backslash function to solve a set of nonlinear equations.
Although this is the most basic nonlinear solver, it is surprisingly powerful. It is based on the NewtonRaphson method in chapter 9.67 of Numerical Recipes in C. In general for well behaved functions and decent initial guesses, its convergence is at least quadratic. However it may fail if the there are local minimums, the condition of the Jacobian is poor or the initial guess is relatively far from the solution. When convergence is negative, it will attempt to backtrack and linesearch. It was validated with fsolve from the MATLAB Optimization Toolbox and IPOPT (https://projects.coinor.org/Ipopt). Please see the help comments and the example.
Note: LSQ curvefit type problems can also be solved using newtonraphson. These are problems where there are many data for a single function, but the coefficients of the function are unknown. Since there is more data than unknowns, and the residual is minimized for each data point, the Jacobian is not square. These problems usually exit with flag 2: "May have converged." and the solution is the best fit in the "leastsquares" sense.
Credit: Moody diagram is from Wikipedia
http://upload.wikimedia.org/wikipedia/commons/8/80/Moody_diagram.jpg
1.6  * allow sparse matrices, replace cond() with condest()


1.5  version 0.4  exit if any dx is nan or inf, allow lsq curvefit type problems. 

1.4  (1) fix need dummy vars dx and convergence to display 0th iteration bug, and (2) if f isnan or isinf need to set lambda2 and f2 before continue bug (3) also reattach moody chart 

1.2  v0.3: Display RCOND and lambda each step. Use ducktyping in funwrapper. Remove Ftyp and F scaling. Use backtracking line search. Output messages, exitflag and min relative step. 

1.1  remove TypicalX and FinDiffRelStep, change feval to evalf & new example solves pipe flow problem using implicit Colebrook equation 
Thomas (view profile)
kang (view profile)
This is a very useful program !!
chenyang (view profile)
RAMI (view profile)
Sir can I get code for three phase newton raphson/backward forward sweep load flow analysis code?
engr.ramzan008@gmail.com
DAO Duc Thu (view profile)
merci bc
Igor Yalovetsky (view profile)
Could you please tell why when b is less tahn zero lambda = slope/(b+sqrt(discriminant))
ali (view profile)
niceone
Naty S (view profile)
Mark Mikofski (view profile)
@Regan Dhakshnamurthy, the latest changes to newtonraphson now allow sparse matrices, and will also quit if your jacobian is singular, but there is no preconditioning. You could hack in to lines: 9899 where it says "% TODO: let user set weights
weight = ones(numel(FUN(x0)),1);" and change the weights to something different, it's not really the same, but this is the matrix used to scale the Jacobian as you can see on line: 112, "Jstar = J./J0; % scale Jacobian" Good luck and let me know if it works and I will make the weights argument callable
Mark Mikofski (view profile)
@ Regan: My guess is no, because if your matrix is so large that inverting is so expensive that you have to resort to pcg instead of just using the backslash operator, then it will be too big for newtonraphson, which also uses the backslash operator. If you want to try it anyway and you have sparse matrices then you'll have to comment out L108 and L181 or change `cond()` to `condest()`. Also you'll have to rewrite your objectives as `Ax  b = 0`, so that the residuals are zero.
Regan Dhakshnamurthy (view profile)
could I substitute this solver to PCG for solving laplacian matrix..
Mark Mikofski (view profile)
There is a minor bug (or feature depending how you look at it). If your Jacobian has any `Inf` or `Nan`, then `cond` throws a confusing exception: "Input to SVD must not contain NaN or Inf." That's because `cond` calls `svd` or "singlevalue decomposition".
I can't remember why I chose `1/cond(J)` instead of `rcond(J)` which would have returned `NaN`. Regardless  this error means your Jacobian has `NaN` or `Inf`, and is therefore badly conditioned.
It's possible that the characteristic `J0` is the culprit b/c this is `1/TYPX` which is a positive nonzero value less than or equal to one (`TYPX = max(abs(x0), 1)`. Therefore if `x0` is **very** large, then `1/x0` will approach zero, causing `Jstar` to contain `Inf`.
EG: `1/(1/(1e999)) > Inf`
That will cause this error  so lessons learned until a patch arrives:
(1) x0 must not be **very** big, EG: x0 < 1e99,
(2) if you get the `svd` error, your Jacobian is badly conditioned, pick better initial guesses and/or better governing equations and
(3) in the future this exception will be caught and return `rcond = NaN`.
Mark Mikofski (view profile)
@Joe The way I've implemented it in the past is a loop around both the fast and slow systems or a nested loop as you described. Either *may* work. EG: system 1 is a chemical kinetics problem, F+O+M > P+M, system 2 is mass transport, dY/dt + u*dY/dx = D*dY2/dx^2, assume transport is slow (diffusion limited) so kinetics reach steady state (SS) instantly  solve kinetics assuming initial concentrations and removing d/dt terms in all equations EG SS solution. Then solve the diffusion problem to find concentrations, now try the kinetics again, does it change the solution? If so repeat until both solutions converge to each other. If problem is super stiff, then they can be completely decoupled, or you might to decrease your timestep to a very small value, to make the steady state assumption true. Note certain systems of equations are inherently unstable, avoid 1st derivative advection problems and backward/forward difference approximations. In general, if your problem is implicitly formulated, it should be stable. Decrease timestep, meshsize, try staggered mesh and always use center difference approximations. Hope that helps! Good luck!
Joe (view profile)
Hi,
This is a very useful program !!
In the comments and ratings as listed below, there was your suggestion using two solution methods to relax the stiffness.
But, I could not make a success in implementing the first solution method, hybrid approach(you have denoted), by iterating a separate calculation for slowest and fastest equations. I have some questions on this approach.
1. Should two equations be solved in its own independent whileloop with different values for tolfun ? This means two whileloops which consists of inner and outer whileloop.
2. If that is the case, is the calculation of solving the fastest equation placed either at inner loop or at outer loop ?
Mark Mikofski (view profile)
The Git repository has been moved from Gist to Github: https://github.com/mikofski/NewtonRaphson. You can report issues there on the issues page. The Gist will be deleted in the near future.
Zoe (view profile)
very useful
Mark Mikofski (view profile)
Fixed a huge bug  the residual function should not be scaled by F(typx), because F(typx) > 0, therefore Ftyp > Inf. Unfortunately previous solution may have been erroneous. Sorry! Latest is in this Gist
https://gist.github.com/mikofski/6381323
and should be online here as soon as TMW reviews it. I've also added backtracking linesearch from Numerical Recipes to improve convergence with poor initial guesses.
Mark Mikofski (view profile)
Ugh, there is an error in funwrapper! Can't read nargout of an anonymous function, but ducktyping works! Just replace
`if abs(nargout(fun))>=2`
with `try` and `else` with `catch`, and voila it works! I will upload a new version when I get a chance, but it is already in the Gist.
https://gist.github.com/mikofski/6381323
Mark Mikofski (view profile)
The pipeflow example uses `IAPWS_IF97()` which can be found on FEX here:
http://www.mathworks.com/matlabcentral/fileexchange/35710iapwsif97functionalformwithnoslip
Mark Mikofski (view profile)
If you have a very old version of MATLAB, you may not be able to use "anonymous functions", but you can replace them with "inline functions".
EG: Anonymous Function
>> f = @(x)cosd(x)
>> f(90)
ans =
0
EG: Inline Function
>> f = inline('cosd(x)','x')
>> f(90)
ans =
0
Mark Mikofski (view profile)
I've uploaded a new version without the optimization only options: FinDiffRelStep is fixed at eps^(1/3) & TypicalX is set to x0, except for zeros which are replaced by ones. Also I improved the example with a pipe flow problem that uses Colebrook and DarcyWeisbach equations, which are implicitly coupled (e.g. nonlinear). Should be online within a week.
For those who can't wait, here is the gist:
https://gist.github.com/mikofski/6381323
Mark Mikofski (view profile)
also evidently typicalX is also not in base MATLAB. :(
Mark Mikofski (view profile)
I just discovered that 'FinDiffRelStop' is not available in base MATLAB; it is only part of the Optimization Toolbox. My apologies. Until I have a chance to update newtonraphson.m please delete 'FinDiffRelStep', eps^(1/3) from oldopts and set dx = eps^(1/3) in the Jacobian function.
Mark Mikofski (view profile)
NOTE: if the Jacobian is badly conditioned, then there will be very flat gradients relative to the steepest gradients, IE: there are very fast changing values and very slow changing values. Another way of saying this is that the problem is "stiff". MATLAB has 3 tools to check the condition of the Jacobian: cond [1], condest [2] & rcond [3]. NOTE: rcond = 1/condest. An even faster estimate of condition is just max(J(:))/min(J(:)). A larger condition means a stiffer system of equations, and possibly harder to solve.
What to do? The classic solution is to use a psuedo (or quasi) steadystate (or equilibrium) assumption. Either (a) solve the fastest equations with everything else held constant first and assume that equations is so fast it jumps straight to that solution, and so it is solved explicitly and can be removed from the problem, or (b) assume the slowest equations are so slow that they're constant, and so can be removed from the problem. If the solutions don't converge take a hybrid approach and iterate between the two extremes: while resnorm>tol & Niter<maxiter, solve(a); solve(b); end.
Another approach is to use damping, in other words limit the size of the change in x from guess to guess. IE: dx_star = min(dx_star, dx_max). A common approach is to scale dx with the norm of the residuals, so that when the norm is big, dx is small and as the problem converges it is damped less and less. IE: dx = dx * tol / resnorm, NOTE: resnorm/tol > +1 and resnorm > tol always.
For an excellent reference on Newton methods see Ch. 23 of Nonlinear Finite Element Methods [4].
[1] http://www.mathworks.com/help/matlab/ref/condest.html
[2] http://www.mathworks.com/help/matlab/ref/cond.html
[3] http://www.mathworks.com/help/matlab/ref/rcond.html
[4] http://www.colorado.edu/engineering/cas/courses.d/NFEM.d/NFEM.Ch23.d/NFEM.Ch23.pdf