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.
}
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}]
}
}