*! version 1.1.2 2/28/95 (STB-24, ssi5.4) program define ridder version 3.1 parse "`*'", parse(" X=") local toload "exp1" while "`1'"!="=" & "`1'"!="returns" & "`1'"!="" { if "`1'"=="X" { if "`toload'"=="exp2" { di in red "X may be specifed only once" exit 198 } local toload "exp2" } else { local `toload' "``toload'' `1'" } mac shift } if "`toload'"!="exp2" { di in red "X must be specified" exit 198 } if "`1'"=="" { error 198 } if "`1'"=="returns" { if "`2'"=="macro" { local type "`3'" local rtnval } else if "`2'"=="exp" { local rtnval "`3'" local type } else error 198 mac shift 3 local progttl "`exp1' X `exp2'" local inline "*" } else { local progttl = substr("`exp2'",2,.) local progttl "`exp1'X`progttl'" local callout "*" } if "`1'" != "=" { di in red "=exp required" exit 198 } mac shift local toload "mval" while "`1'"!="" { if "`1'"=="to" { local toload "mhi" } else if "`1'"=="from" { local toload "mlo" } else if "`1'"=="tol" { local toload "mtol" } else if "`1'"=="ltol" { local toload "mltol" } else { local `toload' "``toload''`1'" } mac shift } if "`mval'"=="" { di in red "= must be followed by expression" exit 198 } if "`mhi'"=="" { di in red "to specifier requried" exit 198 } if "`mlo'"=="" { di in red "from specifier required" exit 198 } tempname flo fhi val lo hi hold ftry ftry2 tol ltol if "`mtol'"=="" & "`mltol'"=="" { scalar `tol' = 1e-6 scalar `ltol' = . } else if "`mltol'"=="" { scalar `tol' = `mtol' scalar `ltol'= . } else if "`mtol'"=="" { scalar `tol' = . scalar `ltol' =`mltol' } else { scalar `tol' = `mtol' scalar `ltol' = `mltol' } scalar `val' = `mval' scalar `lo' = `mlo' scalar `hi' = `mhi' scalar `ltol' = max(`ltol'*abs(`val'),`ltol') local maxiter 64 if `lo'>`hi' { scalar `hold' = `lo' scalar `lo' = `hi' scalar `hi' = `hold' } local ttl : di "Find `progttl'=f() == c=" `val' ", |X_error|<" `tol' /* */ ", |f()-c|<" `ltol' local skip = max(0,int((72-length("`ttl'"))/2)) di in gr _n _skip(`skip') "`ttl'" _n `callout' quietly `exp1' `lo' `exp2' `callout' scalar `flo' = ${`type'} `rtnval' - `val' `inline' scalar `flo' = (`exp1' `lo' `exp2') - `val' `callout' quietly `exp1' `hi' `exp2' `callout' scalar `fhi' = ${`type'} `rtnval' - `val' `inline' scalar `fhi' = (`exp1' `hi' `exp2') - `val' di in gr /* */ _col(17) "lower" _col(39) "upper" _col(59) "weighted" _n /* */ "iteration" _col(17) "bound" _col(27) "f()-c" /* */ _col(39) "bound" _col(49) "f()-c" /* */ _col(59) "midpoint" _col(72) "f()-c" _n /* */ _dup(76) "-" /* verify we bound the solution */ if `flo'*`fhi'>0 | `flo'*`fhi'==. { di in gr %8.0f 0 ". " %9.0g `lo' " " %9.0g `flo' /* */ " " %9.0g `hi' " " %9.0g `fhi' display in red "range does not bound solution" exit 409 } if `flo'==0 { scalar `fhi' = `flo' scalar `hi' = `lo' } else if `fhi'==0 { scalar `flo' = `fhi' scalar `lo' = `hi' } local try = (`lo'+`hi')/2 local lasttry = `try' `callout' quietly `exp1' `try' `exp2' `callout' scalar `ftry' = ${`type'} `rtnval' - `val' `inline' scalar `ftry' = (`exp1' `try' `exp2') - `val' local try = `try'+(`try'-`lo')*sign(`flo'-`fhi')*`ftry' / /* */ sqrt(`ftry'^2 - `flo'*`fhi') `callout' quietly `exp1' `try' `exp2' `callout' scalar `ftry' = ${`type'} `rtnval' - `val' `inline' scalar `ftry' = (`exp1' `try' `exp2') - `val' local iter 0 while `ftry'~=0 &((abs(`lasttry'-`try')>`tol' |abs(`ftry')>`ltol')) { di in gr %8.0f `iter' ". " %9.0g `lo' " " %9.0g `flo' /* */ " " %9.0g `hi' " " %9.0g `fhi' " " /* */ in ye %9.0g `try' " " %9.0g `ftry' local lasttry = `try' if `flo'*`ftry'<0 { scalar `fhi' = `ftry' scalar `hi'=`try' } else { scalar `flo' = `ftry' scalar `lo' = `try' } local try = (`lo'+`hi')/2 `callout' quietly `exp1' `try' `exp2' `callout' scalar `ftry' = ${`type'} `rtnval' - `val' `inline' scalar `ftry' = (`exp1' `try' `exp2') - `val' local try = `try'+(`try'-`lo')*sign(`flo'-`fhi')*`ftry' / /* */ sqrt(`ftry'^2 - `flo'*`fhi') `callout' quietly `exp1' `try' `exp2' `callout' scalar `ftry' = ${`type'} `rtnval' - `val' `inline' scalar `ftry' = (`exp1' `try' `exp2') - `val' local iter = `iter' + 1 if `iter' > `maxiter' { global S_1 = `try' di in red "Convergence not achieved" exit 430 } } di in gr %8.0f `iter' ". " %9.0g `lo' " " %9.0g `flo' /* */ " " %9.0g `hi' " " %9.0g `fhi' " " /* */ in ye %9.0g `try' " " %9.0g `ftry' di in ye _col(57) "|X_error|<=" %9.0g abs(`try'-`lasttry') global S_1 = `try' end exit binsear .002 .005 1e-6 ------------------- Note: specify: rtnval: S_1 to mean $S_1 rtnval: "exp " to mean Logic: $exp will evaluate to Find simul()=f() == c=.22071 1 2 3 4 5 6 7 8 12345678901234567890123456789012345678901234567890123456789012345678901234567890 lower upper iteration bound f()-c bound f()-c midpoint f()-c ---------------------------------------------------------------------------- 123456789 123456789 123456789 123456789 123456789 123456789 123456789 1 2 3 4 5 6 7 8 12345678901234567890123456789012345678901234567890123456789012345678901234567890