-- nderiv.adb compute formulas for numerical derivatives -- order is order of derivative, 1 = first derivative, 2 = second -- points is number of points where value of function is known -- f(x0), f(x1), f(x2) ... f(x points-1) -- term is the point where derivative is computer -- f'(x0) = (1/bh)*( a(0)*f(x0) + a(1)*f(x1) -- + ... + a(points-1)*f(x points-1) -- algorithm: use divided differences to get polynomial p(x) that -- approximates f(x). f(x)=p(x)+error term -- f'(x) = p'(x) + error term' -- substitute xj = x0 + j*h -- substitute x = x0 to get p'(x0) etc with Ada.text_io; use Ada.text_io; with Ada.Numerics; use Ada.Numerics; procedure nderiv is package fltio is new Float_IO(Long_Float); use fltio; type integer_vector is array(integer range <>) of integer; type integer_matrix is array(integer range <>, integer range <>) of integer; norder : integer := 6; npoints : integer := 12; a : integer_vector(0..50); -- coefficients of f(x0), f(x1), ... b : integer; function gcd(a: integer; b : integer) return integer is a1, b1, q, r: integer; begin if a=0 or b=0 then return 1; end if; if abs(a)>abs(b) then a1 := abs(a); b1 := abs(b); else a1 := abs(b); b1 := abs(a); end if; r := 1; while r/=0 loop q := a1 / b1; r := a1 - q * b1; a1 := b1; b1 := r; end loop; return a1; end gcd; procedure deriv(order : integer; npoints : integer; point : integer; a : out integer_vector; bb : out integer) is -- compute the exact coefficients to numerically compute a derivative -- of order 'order' using 'npoints' at ordinate point, -- where order>=1, npoints>=order+1, 0 <= point < npoints, -- 'a' array returned with numerator of coefficients, -- 'bb' returned with denominator of h^order h : integer_vector(0..99); -- coefficients of h x : integer_vector(0..99); -- coefficients of x, variable for differentiating numer : integer_matrix(0..99,0..99); -- numerator of a term numer(term,pos) denom : integer_vector(0..99); -- denominator of coefficient j, k, b : integer; ipower, isum, iat, r : integer; debug : integer := 0; begin if npoints<=order then put_line("ERROR in call to deriv, npoints=" & integer'image(npoints) & "< order=%d+1=" & integer'image(order)); return; end if; for term in 0..npoints-1 loop denom(term) := 1; j := 0; for i in 0..npoints-2 loop if j=term then j := j+1; end if; -- no index j in jth term numer(term,i) := -j; -- from (x-x(j)) denom(term) := denom(term)*(term-j); j := j+1; end loop; end loop; -- end term setting up numerator and denominator if debug>0 then for i in 0..npoints-1 loop put_line("denom(" & integer'image(i) & ")=" & integer'image(denom(i))); for j in 0..npoints-2 loop put_line("numer(" & integer'image(i) & "," & integer'image(j) & ")=" & integer'image(numer(i,j))); end loop; end loop; end if; -- debug -- substitute xj = x0 + j*h, just keep j value in numer(term,i) -- don't need to count x0 because same as x's -- done by above setup -- multiply out polynomial in x with coefficients of h -- starting from (x+j1*h)*(x+j2*h) * ... -- then differentiate d/dx for term in 0..npoints-1 loop x(0) := 1; -- initialize x(1) := 0; h(1) := 0; k := 1; for i in 0..npoints-2 loop for j in 0..k-1 loop -- multiply by (x + j h) h(j) := x(j)*numer(term,i); -- multiply be constant x(j) end loop; for j in 0..k-1 loop -- multiply by x h(j+1) := h(j+1)+x(j); -- multiply by x and add end loop; k := k+1; for j in 0..k-1 loop x(j) := h(j); end loop; x(k) := 0; h(k) := 0; end loop; for j in 0..k-1 loop numer(term,j) := x(j); if debug>0 then put_line("p(x)=" & integer'image(numer(term,j)) & " x^" & integer'image(j) & " at term=" & integer'image(term)); end if; end loop; -- have p(x) for this 'term' -- differentiate 'order' number of times for jorder in 0..order-1 loop for i in 1..k-1 loop -- differentiate p'(x) numer(term,i-1) := i*numer(term,i); end loop; k := k-1; end loop; -- have p^(order) for this term if debug>0 then for i in 0..k-1 loop put_line("p^(" & integer'image(order) & ")(x)=" & integer'image(numer(term,i)) & " x^" & integer'image(i) & " at term=" & integer'image(term)); end loop; end if; end loop; -- end 'term' loop for one 'order', for one 'npoints' -- substitute each point for x, evaluate, get coefficients of f(x(j)) iat := point; -- evaluate at x(iat), substitute for jterm in 0..npoints-1 loop -- get each term ipower := iat; isum := numer(jterm,0); for j in 1..k-1 loop isum := isum + ipower * numer(jterm,j); ipower := ipower * iat; end loop; a(jterm) := isum; if debug>0 then put_line("f^(" & integer'image(order) & ")(x(" & integer'image(iat) & ")) = (1/h^" & integer'image(order) & ") (" & integer'image(a(jterm)) & "/" & integer'image(denom(jterm)) & " f(x(" & integer'image(jterm) & ")) + "); end if; end loop; if debug>0 then new_line; end if; b := 0; for jterm in 0..npoints-1 loop -- clean up each term if a(jterm)/=0 then r := gcd(a(jterm),denom(jterm)); a(jterm) := a(jterm) / r; denom(jterm) := denom(jterm) /r; if denom(jterm)<0 then denom(jterm) := -denom(jterm); a(jterm) := - a(jterm); end if; if denom(jterm)>b then b := denom(jterm); -- largest denominator end if; end if; end loop; for jterm in 0..npoints-1 loop -- check for divisor r := (a(jterm) * b) / denom(jterm); if r * denom(jterm) /= a(jterm) * b then r := gcd(b, denom(jterm)); b := b * (denom(jterm) / r); end if; end loop; for jterm in 0..npoints-1 loop -- adjust for divisor a(jterm) := (a(jterm) * b) / denom(jterm); if debug>0 then put_line("f^(" & integer'image(order) & ")(x(" & integer'image(iat) & "))=(1/" & integer'image(b) & "h^" & integer'image(order) & ")(" & integer'image(a(jterm)) & " f(x("&integer'image(jterm)&")) + "); end if; end loop; -- end computing terms of coefficients bb := b; end deriv; function pow(x : long_float; pwr : integer) return long_float is val : long_float := 1.0; begin for i in 0..pwr-1 loop val := val * x; end loop; return val; end pow; procedure check(order : integer; npoints : integer; point : integer; a : integer_vector; b : integer) is pwr, isum : integer; -- simple test using x^pwr sum, x, term, err, xpoint, deriv : long_float; h : long_float; begin for degree in npoints-2..npoints+1 loop h := 0.125; if order>3 then h := h/2.0; end if; -- cumulative if order>4 then h := h/2.0; end if; if order>5 then h := h/2.0; end if; if order>6 then h := h/2.0; end if; pwr := degree; sum := 0.0; x := 1.0; for i in 0..npoints-1 loop term := long_float(a(i))*pow(x, pwr); sum := sum + term; x := x + h; end loop; sum := sum / long_float(b); sum := sum / pow(h, order); xpoint := 1.0+long_float(point)*h; deriv := 1.0; for i in 0..order-1 loop deriv := deriv * long_float(pwr); pwr := pwr -1; end loop; deriv := deriv * pow(xpoint, pwr); err := deriv - sum; put(" err=" & long_float'image(err) & ", h="); put(h,1,4,0); put(", deriv="); put(deriv,4,1,0); put(", x="); put(xpoint,1,4,0); Put(", pwr=" & integer'image(degree)); new_line; end loop; isum := 0; for i in 0..npoints-1 loop isum := isum + a(I); end loop; if(isum /= 0) then put_Line("ERROR sum a(i) /= 0, is "&integer'image(isum)); end if; end check; begin -- nderiv put_line("Formulas for numerical computation of derivatives "); put_line(" f^(3) means third derivative "); put_line(" f^(2)(x(1)) means second derivative computed at point x(1) "); put_line(" f(x(2)) means function evaluated at point x(2) "); put_line(" Points are x(0), x(1)=x(0)+h, x(2)=x0+2h, ... "); put_line(" The error term gives the power of h and "); put_line(" derivative order with z chosen for maximum value."); put_line(" The error terms, err=, are for test polynomial sum(x^pwr) "); put_line(" with h = 0.125 or smaller, and order = 'pwr' "); for order in 1..norder loop npoints := npoints -2; -- fewer for larger derivatives if npoints<=order+1 then npoints := order + 2; end if; for points in order+1..npoints loop for term in 0..points-1 loop put_line("computing order=" & integer'image(order) & ", npoints=" & integer'image(points) & ", at term=" & integer'image(term)); deriv(order, points, term, a, b); for jterm in 0..points-1 loop if jterm=0 then put_line("f^(" & integer'image(order) & ")(x(" & integer'image(term) & "))=(1/" & integer'image(b) & "h^" & integer'image(order) & ")( " & integer'image(a(jterm)) & " * f(x(" & integer'image(jterm) & ")) + "); elsif jterm=points-1 then put_line(" " & integer'image(a(jterm)) & " * f(x(" & integer'image(jterm) & ") ) + O(h^" & integer'image(points-order) & ")f^(" & integer'image(points) & ")(z) "); else put_line(" " & integer'image(a(jterm)) & " * f(x(" & integer'image(jterm) & ") + "); end if; end loop; -- jterm check(order, points, term, a, b); new_line; end loop; -- term new_line; end loop; -- points new_line; end loop; -- order -- extra checks put_line("order=4, points=9, term=0 "); deriv(4, 9, 0, a, b); check(4, 9, 0, a, b); put_line("order=4, points=9, term=1 "); deriv(4, 9, 1, a, b); check(4, 9, 1, a, b); put_line("order=4, points=10, term=0 "); deriv(4, 10, 0, a, b); check(4, 10, 0, a, b); put_line("order=4, points=10, term=1 "); deriv(4, 10, 1, a, b); check(4, 10, 1, a, b); new_line; put_line("numerical accuracy fails at 4,11 multiple precision required"); new_line; put_line("order=4, points=11, term=0 "); deriv(4, 11, 0, a, b); check(4, 11, 0, a, b); put_line("order=4, points=11, term=1 "); deriv(4, 11, 1, a, b); check(4, 11, 1, a, b); put_line("order=4, points=15, term=0 "); deriv(4, 15, 0, a, b); check(4, 15, 0, a, b); put_line("order=4, points=15, term=1 "); deriv(4, 15, 1, a, b); check(4, 15, 1, a, b); put_line("many formulas checked against Abramowitz and Stegun,"); put_line("Handbook of Mathematical Functions Table 25.2"); end nderiv;