if 0 {
'''Newton-Raphson iterative equation solver'''

This technique can be used to solve many equations of the form f(x)=constant
There is a math procedure called [[''::math::calculus::newtonRaphson func deriv initval'']]
this is bulky in that you need to provide 2 function names to the procedure, and create 2 procs.
Also you need to be able to differentiate the function you are solving.
Here is a simpler code which defines the function in a 'mathematical expression' format
for example [[newtrap pow(x,4) 256]] or [[newtrap x*x*x*x 256]]
will solve  x*x*x*x=256

Newton Raphson uses an iterative scheme to solve the equation, with a better
guess at the solution found from guess = guess - (f(guess)-rhs)/(derivative f(guess))
Newton-Raphson breaks down if a very small derivative is found (eg at x=0 in the x*x*x*x example).
This method uses a numerical differentiation proc (which requires multiple 
evaluations of the function) so is less efficient, but does not require writing 
a proc to evaluate the derivative and will work for any function that can be 
numerically differentiated. NB I have discovered that this technique is called the Secant Method although I think it is morally equivalent to the Newton-Raphson. Mathematicians will immediately produce functions where the actual derivative is not the limit of dy/dx just to upset me of course, but these functions are rarely involved in Physics.
----
 }

 proc derivative {f x} { ;# calculate a numerical derivative for f
	set guess [expr {0.999*$x}]
	set y1 [expr $f]
	set guess [expr {1.001*$x}]
	set y2 [expr $f]
	return [expr {($y2-$y1)/(0.002*$x)}] 
 }

 proc newtrap {f rhs} {
	# convert the lhs of the equation into Newton-Raphson form (lhs-rhs=0)
	# and convert all occurrences of x to the adjustable variable 'guess'
	set f [regsub -all "x" "$f-$rhs" "\$guess"]
	set guess $rhs ;# one possible first guess, not suitable for all equations
	set dd 9999 ;# dd is the change in solution on each iteration; a large value means it has not converged
	set nit 0 ;# limit the number of iterations
	while {$dd>0.00001 || $dd<-0.00001} {
		set dd [expr double([expr $f])/[derivative $f $guess]] ;# change the guess by this
		set guess [expr {$guess-$dd}] ;# new guess
		incr nit
		if {$nit>50	} break ;# time to end
	}
	return $guess
  }

# Examples, solving 3 different equations

  foreach f {"x*x" "x*x-3*x" "pow(x,3)-2*pow(x,2)-5*x" } {
	set x 1
	while {$x<32} {
	    if {[catch {
		set as [newtrap $f $x ]
		puts " solution of $f =$x is at $as"
	    } errmsg]} { puts "at $x $f ($x) $errmsg"}
	    set x [expr {$x+1}]
	}
  }

# more examples solving some inverse trig functions.
  foreach f {asin(x) atan(x) acos(x)} {
  	set x 0.02
  	while {$x<1.000} {
	    if {[catch {
		set as [newtrap $f $x ]
		puts " solution of $f =$x is at $as"
	    } errmsg]} { puts "at $x $f ($x) $errmsg"}
	    set x [expr {$x+.02}]
	}
  }

----
[AM] Interesting idea - but rather than an arbitrary 0.001, it is better to use a factor
that depends on the accuracy - sqrt(epsilon) is a much advocated value - and you need 
to take care of "zero". If the starting estimate is 0, your script will fail.

[GWM] Good points - I was mainly interested in getting solutions to equations involved in map projections such as the Mollweide projection of the earth [http://mathworld.wolfram.com/MollweideProjection.html] which need a solution to the equation:
 2.theta + sin(2*theta) = pi.sin(latitude).

There are a number of other projections requiring similar non-linear equations, and in my case accuracy of 0.001 is fine for an angle in the range 0-pi/2 (latitude never goes above +/-90 deg). I have used this technique successfully for maps of the world, and will wikify the maps when I am happy with them. The projections never (I think) have a zero slope in these equations since zero slope would imply 2 latitudes at the same paper (or canvas) coordinate which would be confusing.

The other errors (derivative calculation at x=0, and some equations require a different choice of starting value to ensure that the equation does not explode in the first few iterations) should be attended to at some point. "Let the User beware!"

I think it would be useful to have a number of mathematical procedures which used a similar 'natural language' interface to defining their inputs. Perhaps a simplified mathematica or maple?

[AM] Don't make me think of the possibilities! :) I have dreamed about something like that for a long time now ...

----
[EKB] Here's a version that does a couple of the things [AM] suggested. It uses a smaller (but user-adjustable) eps and it jitters around the current guess if the derivative is too small. (I also went ahead and called it "secant" even thought I agree -- morally it's Newton-Raphson. But it helps to distinguish from the Newton-Raphson procedure, and let's users know that they don't need to supply a derivative.)

 proc derivative {f x {eps 1.0e-7}} { ;# calculate a numerical derivative for f
        set deltax [expr {$eps * $x}]
        set guess [expr {$x - $deltax}]
	set y1 [expr $f]
	set guess [expr {$x + $deltax}]
	set y2 [expr $f]
	return [expr {($y2-$y1)/(2.0 * $deltax)}]
 }

 proc secant {f rhs {eps 1.0e-7}} {
        set bigeps [expr {sqrt($eps)}]
	# convert the lhs of the equation into Newton-Raphson form (lhs-rhs=0)
	# and convert all occurrences of x to the adjustable variable 'guess'
	set f [regsub -all "x" "$f-$rhs" "\$guess"]
	set guess $rhs ;# one possible first guess, not suitable for all equations
        # Have we made a lucky guess?
        if {abs([expr $f] - $rhs) < $eps} {
            return $guess
        }
	set dd 1.0 ;# dd is the change in solution on each iteration; a value greater than eps means it has not converged
	set nit 0 ;# limit the number of iterations
	while {abs($dd) > $eps} {
               # Jump back and forth, increasing step each time, if we're too close to a small derivative
               set step $bigeps
               while {abs([set deriv [derivative $f $guess]]) < $bigeps * abs([expr $f] - $rhs)} {
                   set guess [expr {$guess + $step}]
                   set step [expr {-2.0 * $step}]
               }
	       set dd [expr double([expr $f])/$deriv] ;# change the guess by this
	       set guess [expr {$guess-$dd}] ;# new guess
	       incr nit
	       if {$nit>50} break ;# time to end
	}
	return $guess
  }

# Problem Examples 
 puts [secant "x * (x-1.0)" 0.5]

 puts [secant "x * x" 0.0]


----
[[[Category Mathematics]|[Category Physics]]]