-- test_lsfit2.adb using stuff from pde_nl33.adb for test with Ada.Text_IO; use Ada.Text_IO; with Ada.Numerics.Long_Elementary_Functions; use Ada.Numerics.Long_Elementary_Functions; with Real_Arrays; -- types real, real_vector, real_matrix use Real_Arrays; -- no higher arrays needed with Integer_Arrays; -- types integer_vector use Integer_Arrays; with Lsfit; use Lsfit; procedure Test_Lsfit2 is Nx : Integer := 15; -- number of data in each dimension Ny : Integer := 15; Xg : Real_Vector(1..Nx); -- x grid, does not need to be uniform Yg : Real_Vector(1..Ny); -- y grid Ug : Real_Matrix(1..Nx,1..Ny); -- computed solution at xyz grid, Ua : Real_matrix(1..Nx,1..Ny); -- known includes boundary Hx, Hy : real; xmin : Real := -1.0; xmax : Real := 1.0; ymin : Real := -1.0; ymax : Real := 1.0; X, Y, U, Ue : Real := 0.0; Last : Boolean := False; Err, Maxerr : Real := 0.0; function Ub_orig(x:Real; y:Real) return real is -- boundary begin return x*x*x*x + 2.0*y*y*y*y + 3.0*x*x*x*y + 4.0*x*y*y*y + 5.0*x*y + 6.0*y*y*y + 7.0*x*x*x + 8.0*x + 9.0*y + 10.0; end Ub_orig; function Ubt(x:Real; y:Real) return real is -- boundary begin return 1.0*x*x*x*x + 2.0*y*y*y*y + 3.0; end Ubt; function ub(x:Real; y:Real) return real is -- data begin return 1.0 + 2.0*x + 3.0*y + 4.0*x*x + 5.0*x*y + 6.0*y*y + 7.0*x*x*x + 8.0*x*x*y + 9.0*x*y*y +10.0*y*y*y + 11.0*x*x*x*x + 12.0*x*x*x*y + 13.0*x*x*y*y + 14.0*x*y*y*y + 15.0*y*y*y*y; end ub; function u1(x:Real; y:Real) return real is -- data begin return 1.0 + 2.0*x + 3.0*y; end u1; function u2(x:Real; y:Real) return real is -- data begin return 1.0 + 2.0*x + 3.0*y + 4.0*x*x + 5.0*x*y + 6.0*y*y; end u2; function u3(x:Real; y:Real) return real is -- data begin return 1.0 + 2.0*x + 3.0*y + 4.0*x*x + 5.0*x*y + 6.0*y*y + 7.0*x*x*x + 8.0*x*x*y + 9.0*x*y*y +10.0*y*y*y; end u3; function u4(x:Real; y:Real) return real is -- data begin return 1.0 + 2.0*x + 3.0*y + 4.0*x*x + 5.0*x*y + 6.0*y*y + 7.0*x*x*x + 8.0*x*x*y + 9.0*x*y*y +10.0*y*y*y + 11.0*x*x*x*x + 12.0*x*x*x*y + 13.0*x*x*y*y + 14.0*x*y*y*y + 15.0*y*y*y*y; end u4; function u5(x:Real; y:Real) return real is -- data begin return 1.0 + 2.0*x + 3.0*y + 4.0*x*x + 5.0*x*y + 6.0*y*y + 7.0*x*x*x + 8.0*x*x*y + 9.0*x*y*y +10.0*y*y*y + 11.0*x*x*x*x + 12.0*x*x*x*y + 13.0*x*x*y*y + 14.0*x*y*y*y + 15.0*y*y*y*y + 16.0*x*x*x*x*x + 17.0*x*x*x*x*y + 18.0*x*x*x*y*y + 19.0*x*x*y*y*y + 20.0*x*y*y*y*y + 21.0*y*y*y*y*y; end u5; function u6(x:Real; y:Real) return real is -- boundary begin return 1.0 + 2.0*x + 3.0*y + 4.0*x*x + 5.0*x*y + 6.0*y*y + 7.0*x*x*x + 8.0*x*x*y + 9.0*x*y*y +10.0*y*y*y + 11.0*x*x*x*x + 12.0*x*x*x*y + 13.0*x*x*y*y + 14.0*x*y*y*y + 15.0*y*y*y*y + 16.0*x*x*x*x*x + 17.0*x*x*x*x*y + 18.0*x*x*x*y*y + 19.0*x*x*y*y*y + 20.0*x*y*y*y*y + 21.0*y*y*y*y*y + 22.0*x*x*x*x*x*x + 23.0*x*x*x*x*x*y + 24.0*x*x*x*x*y*y + 25.0*x*x*x*y*y*y + 26.0*x*x*y*y*y*y + 27.0*x*y*y*y*y*y + 28.0*y*y*y*y*y*y; end u6; package Real_IO is new Float_IO(Real); use Real_IO; begin Put_Line("test_lsfit2.adb"); hx := (xmax-xmin)/real(nx-1); Put_Line("nx="&Integer'Image(nx)&", step="&Real'Image(hx)); for i in 1..nx loop xg(i) := xmin+Real(i-1)*hx; end loop; hy := (ymax-ymin)/real(ny-1); Put_Line("ny="&Integer'Image(ny)&", step="&Real'Image(hy)); for ii in 1..ny loop yg(ii) := ymin+Real(ii-1)*hy; end loop; for K in 1..6 loop Fit_Init(K, 2); -- Kth power, two independent variables Last := False; for i in 1..nx loop X := Xg(i); for ii in 1..ny loop Y := Yg(ii); if K=1 then U := U1(X,Y); end if; if K=2 then U := U2(X,Y); end if; if K=3 then U := U3(X,Y); end if; if K=4 then U := U4(X,Y); end if; if K=5 then U := U5(X,Y); end if; if K=6 then U := U6(X,Y); end if; if i=nx and ii=ny then Last := True; end if; Fit2_Data(X, Y, U, Last); end loop; end loop; -- built, now check Maxerr := 0.0; for i in 1..nx loop X := Xg(I); for ii in 1..ny loop Y := Yg(ii); if K=1 then U := U1(X,Y); end if; if K=2 then U := U2(X,Y); end if; if K=3 then U := U3(X,Y); end if; if K=4 then U := U4(X,Y); end if; if K=5 then U := U5(X,Y); end if; if K=6 then U := U6(X,Y); end if; Ue := Eval2_poly(X, Y); Err := abs(U-Ue); if Err>Maxerr then Maxerr := Err; end if; end loop; end loop; Put_Line("power "&Integer'Image(K)&" test, maxerr="&Real'Image(Maxerr)); end loop; Put_Line("just boundary to data, power 4"); Fit_Init(4, 2); -- 4th power, 2 dimension Last := False; for i in 1..nx loop X := xg(i); for ii in 1..ny loop Y := yg(ii); Ua(i,ii) := u4(X,Y); -- to check answers inside boundary if i=1 or i=nx or ii=1 or ii=ny then Ug(i,ii) := u4(X,Y); -- just give boundary to fit if i=nx and ii=ny then Last := True; end if; Fit2_Data(X, Y, Ug(i,ii), Last); else Ug(i,ii) := 0.0; -- set up for next test end if; end loop; end loop; New_Line; -- check Maxerr := 0.0; for i in 1..nx loop X := xg(i); for ii in 1..ny loop Y := yg(ii); if i=1 or i=nx or ii=1 or ii=ny then U := u4(X,Y); Ue := Eval2_poly(X, Y); Err := abs(U-Ue); if Err>Maxerr then Maxerr := Err; end if; end if; end loop; end loop; Put_Line("just boundary test4, maxerr="&Real'Image(Maxerr)); New_Line; -- check Maxerr := 0.0; for i in 2..Nx-1 loop X := xg(i); for ii in 2..Ny-1 loop Y := yg(ii); U := u4(X,Y); Ue := Eval2_poly(X, Y); Err := abs(U-Ue); if Err>Maxerr then Maxerr := Err; end if; end loop; end loop; Put_Line("all internal test4, maxerr="&Real'Image(Maxerr)); New_Line; Maxerr := 0.0; Put_Line("check Bpoly power 4, dimension 2"); Fit_Init(4,2); -- 4 power, 2 dimension Fit2_Bpoly(Nx, Ny, Xg, Yg, Ug); -- only boundary of Ug used -- test fit on boundary Maxerr := 0.0; for i in 1..nx loop X := Xg(i); for ii in 1..ny loop Y := Yg(ii); if i=1 or i=nx or ii=1 or ii=ny then U := Eval2_Poly(X, Y); Err := abs(U4(X, Y)-U); if Err>Maxerr then Maxerr := Err; end if; end if; end loop; end loop; Put_Line("Bpoly fit, maxerr="&Real'Image(Maxerr)); New_Line; -- test fit internal Maxerr := 0.0; for i in 2..Nx-1 loop X := Xg(i); for ii in 2..Ny-1 loop Y := Yg(ii); U := Eval2_Poly(X, Y); Err := abs(U4(X, Y)-U); if Err>Maxerr then Maxerr := Err; end if; end loop; end loop; Put_Line("internal Bpoly fit, maxerr="&Real'Image(Maxerr)); New_Line; Put_Line("test_lsfit2.adb finished"); end Test_Lsfit2;