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 [L1 ] 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 ...