if d3=0: .01→d3
(√(abs(expr(ffun)/(d3/(h1^2)))))/10→hh
if hh=0: .01→hh
endif
clrio
disp "Starting interval hh: "&string(hh)
newmat(ntab,ntab)→amat
(expr(fphh)-expr(fmhh))/(2*hh)→amat[1,1]
big→err
for i,2,ntab
hh/con→hh
(expr(fphh)-expr(fmhh))/(2*hh)→amat[1,i]
con2→fac
for j,2,i
(amat[j-1,i]*fac-amat[j-1,i-1 ])/(fac-1)→amat[j,i]
con2*fac→fac
max(abs(amat[j,i]-amat[j-1,i]),abs(amat[j,i]-amat[j-1,i-1]))→errt
if errt≤err then
errt→err
amat[j,i]→dest
endif
endfor
if abs(amat[i,i]-amat[i-1,i-1])≥ safe*err then
disp "SAFE limit exit"
exit
endif
endfor
disp "dy/dx estimate: "&string(dest)
disp "error: "&string(err)
disp "amat[] column: "&string(i)
disp "fac: "&string(fac)
pause
disphome
endprgm
General comments on numerical differentiation
Any numerical derivative routine is going to suffer accuracy problems where the function is changing
rapidly. This includes asymptotes. For example, suppose that we try to find the derivative of
tan(1.5707). This is very close to the asymptote of pi/2 = 1.57079632679. nder2() returns an answer of
-1.009...E6, which is clearly wrong. The problem is that the starting interval brackets the asypmtote.
nder1() returns an answer of -1.038..E7, but at least also returns an error of 3.427E8, so we know
something is wrong.
The function must be continuous on the sample interval for both nder1() and nder2().
6 - 34
To Next Page
Loading...
