-- nuderiv_test_exp.adb non uniformly spaced derivative coefficients with Ada.Text_IO; use Ada.Text_IO; with Ada.Numerics; use Ada.Numerics; with Ada.Numerics.Long_Elementary_Functions; use Ada.Numerics.Long_Elementary_Functions; with real_Arrays; use real_Arrays; with nuderiv; -- procedure under test with rderiv; -- for comparison procedure nuderiv_test_exp is npoint : integer; n : Integer := 12; X : real_vector(1..n); Y : real_vector(1..n); Yp : real_vector(1..n); Ypp : real_vector(1..n); Yppp : real_vector(1..n); Ypppp : real_vector(1..n); Yppppp : real_vector(1..n); C, Cu : real_vector(1..n); nuyp, nyp : real; -- any order Xmin, Xmax, hx : real; -- for comparison test Err1, Err2, Err3, Err4 , Err5 : Real; Max_Err1, Max_Err2, Max_Err3, Max_Err4, Max_Err5 : Real := 0.0; nx : integer; h : Real; package Fltio is new Float_Io(Long_Float); use Fltio; package Intio is new integer_Io(integer); use Intio; -- test function and check functions function yf5(x : real) return real is -- test function 5 begin return exp(x); end yf5; function ypf5(x : real) return real is -- derivative of test function 5 begin return exp(x); end ypf5; function yppf5(x : real) return real is -- 2nd derivative of test function 5 begin return exp(x); end yppf5; function ypppf5(x : real) return real is -- 3rd derivative of test function5 begin return exp(x); end ypppf5; function yppppf5(x : real) return real is -- 4th derivative of test function5 begin return exp(x); end yppppf5; function ypppppf5(x : real) return real is -- 5th derivative of test function5 begin return exp(x); end ypppppf5; begin put_line("nuderiv_test_exp.adb running "); npoint := n; H := 0.1; for I in 1..N loop X(I) := H; H := H*1.33; -- increasing spacing end loop; put_line(" testing npoint="&integer'Image(npoint)); for i in 1..npoint loop Y(i) := yf5(x(i)); Yp(i) := ypf5(x(i)); Ypp(i) := yppf5(x(i)); Yppp(i) := ypppf5(x(i)); Ypppp(i) := yppppf5(x(i)); Yppppp(i) := ypppppf5(x(i)); put("x("); put(I,1); put(")="); put(X(i),3,3); new_line; Put(" y="); Put(Y(i),3,3); Put(", y'="); Put(Yp(i),3,3); Put(", y''="); Put(Ypp(i),3,3); Put(", y'''="); Put(Yppp(i),3,3); Put(", y''''="); Put(Ypppp(i),3,3); Put(", y'''''="); Put(Yppppp(i),3,3); new_line; end loop; for order in 1..npoint-1 loop for point in 1..npoint loop Put_line("testing order="&integer'image(order)& ", npoint="&integer'image(npoint)& ", point="&integer'image(point)); nuderiv(order, npoint, point, X, C); for i in 1..npoint loop put("nu c("); put(i,1); put(")="); put(C(i),3,3); new_line; end loop; nuyp := 0.0; for i in 1..npoint loop nuyp := nuyp + C(i)*Y(i); end loop; if order=1 then put_line("order="&integer'image(order)& ", npoint="&integer'Image(npoint)& ", point="&Integer'Image(point)); put(" yp="); put(yp(point),3,3); put(", nuyp="); put(nuyp,3,3); put(", err1="); err1 := yp(point)-nuyp; put(err1,3,3); new_line; if abs(err1) > Max_Err1 then Max_Err1 := abs(err1); end if; elsif order=2 then put_line("order="&integer'image(order)& ", npoint="&integer'Image(npoint)& ", point="&Integer'Image(point)); put(" ypp="); put(ypp(point),3,3); put(", nuyp="); put(nuyp,3,3); put(", err2="); err2 := ypp(point)-nuyp; put(err2,3,3); new_line; if abs(err2) > Max_Err2 then Max_Err2 := abs(err2); end if; put(err2,3,3); new_line; elsif order=3 then put_line("order="&integer'image(order)& ", npoint="&integer'Image(npoint)& ", point="&Integer'Image(point)); put(" yppp="); put(yppp(point),3,3); put(", nuyp="); put(nuyp,3,3); put(", err3="); err3 := yppp(point)-nuyp; put(err3,3,3); new_line; if abs(err3) > Max_Err3 then Max_Err3 := abs(err3); end if; put(err3,3,3); new_line; elsif order=4 then put_line("order="&integer'image(order)& ", npoint="&integer'Image(npoint)& ", point="&Integer'Image(point)); put(" ypppp="); put(ypppp(point),3,3); put(", nuyp="); put(nuyp,3,3); put(", err4="); err4 := ypppp(point)-nuyp; put(err4,3,3); new_line; if abs(err4) > Max_Err4 then Max_Err4 := abs(err4); end if; put(err4,3,3); new_line; elsif order=5 then put_line("order="&integer'image(order)& ", npoint="&integer'Image(npoint)& ", point="&Integer'Image(point)); put(" yppppp="); put(yppppp(point),3,3); put(", nuyp="); put(nuyp,3,3); put(", err5="); err5 := yppppp(point)-nuyp; put(err5,3,3); new_line; if abs(err5) > Max_Err5 then Max_Err5 := abs(err5); end if; put(err5,3,3); new_line; end if; end loop; end loop; New_Line; Put_Line("Compare accuracy of nuderiv and nderiv on uniform spaced"); Xmin := 0.0; Xmax := 1.0; nx := 5; Hx := (Xmax-Xmin)/real(nx-1); for i in 1..nx loop X(i) := Xmin + Hx*Real(i-1); Y(i) := X(i)**4 + 2.0*X(I)**3 + 3.0*X(I)**2 + 4.0*X(I) + 5.0; Put_Line("x="&Real'Image(X(i))&", y="&Real'Image(Y(i))); end loop; for i in 1..nx loop -- point nuderiv(4, nx, i, X, C); rderiv( 4, nx, i, Hx, Cu); nuyp := 0.0; nyp := 0.0; Put_Line("testing order 4, 5 points, at point i="&Integer'Image(i)); for j in 1..nx loop nuyp := nuyp + Cu(j)*Y(j); nyp := nyp + C(j)*Y(j); end loop; -- j Put("nu="); Put(nuyp,2,3,0); Put(", nd="); Put(nyp,2,3,0); Put(", actual="); Put(24.0,2,3,0); Put_Line(", diff="&Real'Image(Nuyp-Nyp)); end loop; -- i Put_Line("first order max_err1="&Real'Image(Max_Err1)); Put_Line("second order max_err2="&Real'Image(Max_Err2)); Put_Line("third order max_err3="&Real'Image(Max_Err3)); Put_Line("fourth order max_err4="&Real'Image(Max_Err4)); Put_Line("fifth order max_err5="&Real'Image(Max_Err5)); end nuderiv_test_exp;