-- pde_nl33.adb non linear PDE, second order, second degree, two dimension -- with many comments and much checking -- did not have to be uniform grid -- could have used sparse matrix -- this is a teaching example -- -- The PDE to be solved is: -- -- u(x,y,z)*uxxx(x,y,z)+2*uyyy(x,y,z)*uzzz(x,y,z)+ -- 3*ux(x,y,z)*uy(x,y,z)*uz(x,y,z) = f(x,y,z) -- -- f(x,y,z) = 24*(x^4+2*y^4+3*z^4+4*x*y*z+5*x*y+z+2)*x+6912*y*z+ -- 3*(4*x^3+4*y*z+5*y)*(8*y^3+4*x*z+5*x)*(12*z^3+4*x*y+1) -- -- Boundary conditions computed from -- ub(x,y,z) = x^4+2*y^4+3*z^4+4*x*y*z+5*x*y+z+2 -- On cube -1 To 1 in three dimensions with Ada.Text_IO; use Ada.Text_IO; with Ada.Numerics.Long_Elementary_Functions; use Ada.Numerics.Long_Elementary_Functions; with Ada.Calendar; use Ada.Calendar; with Real_Arrays; -- types real, real_vector, real_matrix use Real_Arrays; with Integer_Arrays; -- types integer_vector use Integer_Arrays; with deriv; with rderiv; with Simeq_Newton5; with Array3d; -- type real_matrix3d use Array3d; with Lsfit; -- least square fit for initial guess use Lsfit; procedure pde_nl33 is nx : Integer := 5; -- number of coefficients for discrete derivative ny : Integer := 5; nz : Integer := 5; dof : Integer := (nx-2)*(ny-2)*(nz-2); -- number of equations and -- number of linear unknowns nlin : Integer := dof*Dof*dof; -- number of possible nonlinear terms, unlin : Integer := 0; -- actual nonlinear terms used ntot : Integer := dof+nlin; -- total number of terms, includes nonlinear A : Real_Matrix(1..dof,1..ntot); -- nonlinear system of equations Yy : Real_Vector(1..dof); -- RHS of equations Xs : Real_Vector(1..ntot); -- solution being computed, not include bound xg : Real_Vector(1..nx); -- x grid, does not need to be uniform yg : Real_Vector(1..ny); -- y grid zg : Real_Vector(1..nz); -- z grid ug : Real_Matrix3d(1..nx,1..ny,1..nz); -- computed solution at xyz grid, -- boundary needed var1 : Integer_Vector(1..ntot); -- for nonlinear simeq var2 : Integer_Vector(1..ntot); -- for nonlinear simeq var3 : Integer_Vector(1..ntot); -- for nonlinear simeq vari1 : Integer_Vector(1..ntot); -- for nonlinear simeq vari2 : Integer_Vector(1..ntot); -- for nonlinear simeq Ua : Real_matrix3d(1..nx,1..ny,1..nz); -- known includes boundary cx : Real_Matrix(1..nx,1..nx); -- discrete first derivative cxxx : Real_Matrix(1..nx,1..nx); -- discrete third derivative hx : real; cy : Real_Matrix(1..ny,1..ny); -- discrete first derivative cyyy : Real_Matrix(1..ny,1..ny); -- discrete third derivative hy : real; cz : Real_Matrix(1..nz,1..nz); -- discrete first derivative czzz : Real_Matrix(1..nz,1..nz); -- discrete third derivative hz : real; err, avgerr, Maxerr, x, y, z : real; eps : real := 1.0e-4; xmin : Real := -1.0; xmax : Real := 1.0; ymin : Real := -1.0; ymax : Real := 1.0; zmin : Real := -1.0; zmax : Real := 1.0; Xnames : array(1..dof) of String(1..4); -- nx<=9, ny<=9, nz<=9 -- e.g. ("X111" ... "X333"); Start_Time, End_Time : Duration; Debug : Boolean := False; function S(i:Integer; ii:Integer; iii:Integer) return Integer is begin return (i-2)*(ny-2)*(nz-2)+(ii-2)*(nz-2)+(iii-2)+1; --make single subscript end S; -- definition of PDE to be solved function f(x:Real; y:Real; z:Real) return real is -- forcing function begin return 24.0*(x*x*x*x+2.0*y*y*y*y+3.0*z*z*z*z+4.0*x*y*z+ 5.0*x*y+z+2.0)*x+6912.0*y*z+ 3.0*(4.0*x*x*x+4.0*y*z+5.0*y)* (8.0*y*y*y+4.0*x*z+5.0*x)*(12.0*z*z*z+4.0*x*y+1.0); end f; function ub(x:Real; y:Real; z:Real) return real is -- boundary begin return x*x*x*x+2.0*y*y*y*y+3.0*z*z*z*z+4.0*x*y*z+5.0*x*y+z+2.0; end ub; package Real_IO is new Float_IO(Real); use Real_IO; package intio is new integer_io(integer); use intio; procedure Build_Vars is -- dof, then: find_term based on PDE begin for i in 1..dof loop var1(i) := i; -- linear terms, required var2(i) := 0; var3(i) := 0; vari1(i) := 0; vari2(i) := 0; end loop; -- unlin grows dof+1..unlin end Build_Vars; function Find_Term(T1, T2, T3 : Integer) return Integer is -- term index index : Integer := 1; V1 : Integer := T1; V2 : Integer := T2; V3 : Integer := T3; V : Integer; begin -- must have largest first for uniqueness if V2>V1 then V:=V1; V1:=V2; V2:=V; end if; if V3>V1 then V:=V1; V1:=V3; V3:=V; end if; if V3>V2 then V:=V2; V2:=V3; V3:=V; end if; while V1/=Var1(index) or V2/=Var2(index) or V3/=Var3(index) loop if index>dof+unlin then Var1(index) := V1; Var2(index) := V2; Var3(index) := V3; unlin := unlin+1; exit; end if; index := index+1; end loop; return index; end Find_Term; procedure Compute_Discrete_Derivatives is T : Real_Vector(1..nx+ny); begin -- first subscript is location i for cx; ii for cy; iii for cz -- second subscript is point j for cx; jj for cy; jjj for cz for i in 1..nx loop rderiv(1, nx, i, hx, T); -- may be non uniform grid, nuderiv for j in 1..nx loop cx(i,j) := T(j); -- discrete derivative coefficients if Debug then Put_Line("cx("&Integer'Image(i)&","&Integer'Image(j)&")="& Real'Image(cx(i,j))); end if; -- debug end loop; rderiv(3, nx, i, hx, T); for j in 1..nx loop cxxx(i,j) := T(j); if Debug then Put_Line("cxxx("&Integer'Image(i)&","&Integer'Image(j)&")="& Real'Image(cxxx(i,j))); end if; -- debug end loop; end loop; for ii in 1..ny loop rderiv(1, ny, ii, hy, T); -- may be non uniform grid, nuderiv for jj in 1..ny loop cy(ii,jj) := T(jj); if Debug then Put_Line("cy("&Integer'Image(ii)&","&Integer'Image(jj)&")="& Real'Image(cy(ii,jj))); end if; -- debug end loop; rderiv(3, ny, ii, hy, T); for jj in 1..ny loop cyyy(ii,jj) := T(jj); if Debug then Put_Line("cyyy("&Integer'Image(ii)&","&Integer'Image(jj)&")="& Real'Image(cyyy(ii,jj))); end if; -- debug end loop; end loop; for iii in 1..nz loop rderiv(1, nz, iii, hz, T); -- may be non uniform grid, nuderiv for jjj in 1..nz loop cz(iii,jjj) := T(jjj); if Debug then Put_Line("cz("&Integer'Image(iii)&","&Integer'Image(jjj)&")="& Real'Image(cz(iii,jjj))); end if; -- debug end loop; rderiv(3, nz, iii, hz, T); for jjj in 1..nz loop czzz(iii,jjj) := T(jjj); if Debug then Put_Line("czzz("&Integer'Image(iii)&","&Integer'Image(jjj)&")="& Real'Image(czzz(iii,jjj))); end if; -- debug end loop; end loop; end Compute_Discrete_Derivatives; procedure Build_Equations is ia, k : Integer; t : Real; -- coefficient of k term being computed begin -- initialize A and Yy for i in 1..dof loop for j in 1..ntot loop A(i,j) := 0.0; end loop; Yy(i) := 0.0; end loop; -- compute entries in nonlinear equation matrix A*Xs=Yy -- method -- u(x,y,z)*uxxx(x,y,z)+2*uyyy(x,y,z)*uzzz(x,y,z)+ -- 3*ux(x,y,z)*uy(x,y,z)*uz(x,y,z) = f(x,y,z) -- becomes system of non linear equations Put_line("compute non linear A matrix"); for i in 2..nx-1 loop -- xg inside boundary x := xg(i); for ii in 2..ny-1 loop -- yg inside boundary y := yg(ii); for iii in 2..nz-1 loop -- zg inside boundary z := zg(iii); ia := S(i,Ii,iii); -- equation index Yy(ia) := f(x,Y,z); -- RHS -- term u(x,y,z)*uxxx(x,y,z) x at i, y at ii, z at iii for j in 1..nx loop -- ug(j,ii,iii) cx(i,j) if j=1 or j=nx then -- boundary in x direction k := Find_Term(S(i,ii,iii),0,0); t := (cxxx(i,j)*ug(j,ii,iii)); A(ia,k) := A(ia,k) + t; else -- j k := Find_Term(S(i,ii,iii),S(j,ii,iii),0); t := cxxx(i,j); A(ia,k) := A(ia,k) + t; end if; -- j end loop; -- j -- term 2*uyyy(x,y,z)*uzzz(x,y,z) for jj in 1..ny loop -- ug(i,jj,iii) cy(ii,jj) if jj=1 or jj=ny then -- boundary in y direction for jjj in 1..nz loop -- ug(i,ii,jjj) cz(iii,jjj) if jjj=1 or jjj=nz then -- boundary in z direction k := 0; t := 2.0*(cyyy(ii,jj)*ug(i,jj,iii))* (czzz(iii,jjj)*ug(i,ii,jjj)); Yy(ia) := Yy(ia) - t; else -- jjj k := Find_Term(S(i,ii,jjj),0,0); t := 2.0*czzz(iii,jjj)*(cyyy(ii,jj)*ug(i,jj,iii)); A(ia,k) := A(ia,k) + t; end if; -- jjj end loop; -- jjj else -- jj for jjj in 1..nz loop -- ug(i,ii,jjj) cz(iii,jjj) if jjj=1 or jjj=nz then -- boundary in z direction k := Find_Term(S(i,jj,iii),0,0); t := 2.0*cyyy(ii,jj)*(czzz(iii,jjj)*ug(i,ii,jjj)); A(ia,k) := A(ia,k) + t; else -- jjj k := Find_Term(S(i,jj,iii),S(i,ii,jjj),0); t := 2.0*cyyy(ii,jj)*czzz(iii,jjj); A(ia,k) := A(ia,k) + t; end if; -- jjj end loop; -- jjj end if; -- jj end loop; -- jj -- term 3*ux(x,y,z)*uy(x,y,z)*uz(x,y,z) for j in 1..nx loop -- ug(j,ii,iii) cx(i,j) if j=1 or j=nx then -- boundary in x direction for jj in 1..ny loop -- ug(i,jj,iii) cy(ii,jj) if jj=1 or jj=ny then -- boundary in y direction for jjj in 1..nz loop -- ug(i,ii,jjj) cz(iii,jjj) if jjj=1 or jjj=nz then -- boundary in z direction k := 0; t := 3.0*(cx(i,j)*ug(j,ii,iii))* (cy(ii,jj)*ug(i,jj,iii))* (cz(iii,jjj)*ug(i,ii,jjj)); Yy(ia) := Yy(ia) - t; else -- jjj k := Find_Term(S(i,ii,jjj),0,0); t := 3.0*cz(iii,jjj)*(cx(i,j)*ug(j,ii,iii))* (cy(ii,jj)*ug(i,jj,iii)); A(ia,k) := A(ia,k) + t; end if; -- jjj end loop; -- jjj else -- jj for jjj in 1..nz loop -- ug(i,ii,jjj) cz(iii,jjj) if jjj=1 or jjj=nz then -- boundary in z direction k := Find_Term(S(i,jj,iii),0,0); t := 3.0*cy(ii,jj)*(cx(i,j)*ug(j,ii,iii))* (cz(iii,jjj)*ug(i,ii,jjj)); A(ia,k) := A(ia,k) + t; else -- jjj k := Find_Term(S(i,jj,iii),S(i,ii,jjj),0); t := 3.0*cy(ii,jj)*cz(iii,jjj)*(cx(i,j)*ug(j,ii,iii)); A(ia,k) := A(ia,k) + t; end if; -- jjj end loop; -- jjj end if; -- jj end loop; -- jj else -- j for jj in 1..ny loop -- ug(i,jj,iii) cy(ii,jj) if jj=1 or jj=ny then -- boundary in y direction for jjj in 1..nz loop -- ug(i,ii,jjj) cz(iii,jjj) if jjj=1 or jjj=nz then -- boundary in z direction k := Find_Term(S(j,ii,iii),0,0); t := 3.0*cx(i,j)*(cy(ii,jj)*ug(i,jj,iii))* (cz(iii,jjj)*ug(i,ii,jjj)); A(ia,k) := A(ia,k) + t; else -- jjj k := Find_Term(S(j,ii,iii),S(i,ii,jjj),0); t := 3.0*cx(i,j)*cz(iii,jjj)*(cy(ii,jj)*ug(i,jj,iii)); A(ia,k) := A(ia,k) + t; end if; -- jjj end loop; -- jjj else -- jj for jjj in 1..nz loop -- ug(i,ii,jjj) cz(iii,jjj) if jjj=1 or jjj=nz then -- boundary in z direction k := Find_Term(S(j,ii,iii),S(i,jj,iii),0); t := 3.0*cx(i,j)*cy(ii,jj)*(cz(iii,jjj)*ug(i,ii,jjj)); A(ia,k) := A(ia,k) + t; else -- jjj k := Find_Term(S(j,ii,iii),S(i,jj,iii),S(i,ii,jjj)); t := 3.0*cx(i,j)*cy(ii,jj)*cz(iii,jjj); A(ia,k) := A(ia,k) + t; end if; -- jjj end loop; -- jjj end if; -- jj end loop; -- jj end if; -- j end loop; --j end loop; -- iii end loop; -- ii end loop; -- i end Build_Equations; -- finished A and Y for all dof equations procedure Print_AY is begin Put_Line("nonlinear matrix A, zero entries not printed"); for i in 1..dof loop for j in 1..dof+unlin loop if A(i,j) /= 0.0 then Put_Line("i="&Integer'Image(i)&", j="&Integer'Image(j)&", A(i,j)="& Real'Image(A(i,j))); end if; end loop; end loop; New_Line; Put_Line("Y computed RHS "); for i in 1..dof loop Put_Line("Y("&Integer'Image(i)&")="&Real'Image(Yy(i))); end loop; New_Line; end Print_AY; procedure check_soln(U : Real_Matrix3d) is -- check when solution unknown smaxerr, err, x, y, z, uv, ux, uxxx, uy, uyyy, uz, uzzz : Real; begin Put_line("check_soln"); smaxerr := 0.0; for i in 2..nx-1 loop -- through dof equations x := xg(I); for ii in 2..ny-1 loop y := yg(ii); for iii in 2..nz-1 loop z := zg(iii); err := -f(x,y,z); -- add other side of equation terms uv := U(i,ii,iii); -- solution being checked ux := 0.0; uxxx := 0.0; for j in 1..nx loop -- compute numerical derivative at (i,ii,iii) ux := ux + cx(i,j)*U(j,ii,iii); uxxx := uxxx + cxxx(i,j)*U(j,ii,iii); end loop; -- j uy := 0.0; uyyy := 0.0; for jj in 1..ny loop uy := uy + cy(ii,jj)*U(i,jj,iii); uyyy := uyyy + cyyy(ii,jj)*U(i,jj,iii); end loop; -- jj uz := 0.0; uzzz := 0.0; for jjj in 1..nz loop uz := uz + cz(iii,jjj)*U(i,ii,jjj); uzzz := uzzz + czzz(iii,jjj)*U(i,ii,jjj); end loop; -- jj err := err + uv*uxxx; err := err + 2.0*uyyy*uzzz; err := err + 3.0*ux*uy*uz; if abs(err) > smaxerr then smaxerr := abs(err); end if; end loop; -- iii Z end loop; -- ii Y end loop; -- i X Put_line("check_soln max error="&Real'Image(smaxerr)); end Check_Soln; procedure Print_Equations(n : integer; unlin : integer; var1 : integer_vector; var2 : integer_vector; var3 : integer_vector; vari1 : integer_vector; vari2 : Integer_Vector) is begin Put_Line("system of equations to be solved, i=1.."&Integer'Image(n)); for j in 1..n loop Put_line("A(i,"&Integer'Image(j)&")*"&Xnames(j)&"+"); end loop; for j in n+1..n+unlin loop Put("A(i,"&Integer'Image(j)&")"); if Var1(j)>0 then Put("*"&Xnames(Var1(j))); end if; if Var2(j)>0 then Put("*"&Xnames(Var2(j))); end if; if Var3(j)>0 then Put("*"&Xnames(Var3(j))); end if; if Vari1(j)>0 then Put("/("&Xnames(Vari1(j))); end if; if Vari2(j)>0 then Put("*"&Xnames(Vari2(j))); end if; if Vari1(j)>0 or Vari2(j)>0 then Put(")"); end if; if J /= n+unlin then Put_Line("+"); else New_Line; end if; end loop; Put_Line(" = Y(i)"); New_Line; end Print_Equations; procedure Check_Equations(n : integer; unlin : integer; A : Real_Matrix; Y : Real_Vector; var1 : integer_vector; var2 : integer_vector; var3 : integer_vector; vari1 : integer_vector; vari2 : Integer_Vector; X : in out Real_vector) is resid : Real := 0.0; begin for i in 1..n loop Put_Line(Xnames(i)&" = "&Real'Image(X(i))); end loop; for j in n+1..n+unlin loop -- set non linear X X(j) := 1.0; if var1(j) > 0 then X(j) := X(j)*X(var1(j)); end if; if var2(j) > 0 then X(j) := X(j)*X(var2(j)); end if; if var3(j) > 0 then X(j) := X(j)*X(var3(j)); end if; if vari1(j) > 0 then X(j) := X(j) / X(vari1(j)); end if; if vari2(j) > 0 then X(j) := X(j) / X(vari2(j)); end if; end loop; -- compute residual, abs sum of all equations for i in 1..n loop for j in 1..n+unlin loop resid := resid + A(i,j)*X(j); end loop; -- j resid := resid - Yy(i); resid := abs(resid); end loop; -- i Put_Line("Check_Equation residual = "&Real'Image(resid)); New_Line; end Check_Equations; procedure Make_Xnames is begin for i in 2..nx-1 loop for ii in 2..ny-1 loop for iii in 2..nz-1 loop Xnames(S(i,ii,iii)) := "X"&Integer'Image(i-1)(2..2)& Integer'Image(ii-1)(2..2)& Integer'Image(iii-1)(2..2); Put_Line(Xnames(S(i,ii,iii))); end loop; end loop; end loop; end Make_Xnames; begin Put_Line("pde_nl33.adb running "); Put_Line("nonlinear second order, second degree, two dimension"); Put_Line("The PDE to be solved for u(x) is:"); Put_Line(" "); Put_Line("u(x,y,z)*uxxx(x,y,z)+2*uyyy(x,y,z)*uzzz(x,y,z)+"); Put_Line("3*ux(x,y,z)*uy(x,y,z)*uz(x,y,z) = f(x,y,z)"); Put_Line(" "); Put_Line("f(x,y,z)=24*(x^4+2*y^4+3*z^4+4*x*y*z+5*x*y+z+2)*x+6912*y*z+"); Put_Line(" 3*(4*x^3+4*y*z+5*y)*(8*y^3+4*x*z+5*x)*(12*z^3+4*x*y+1)"); Put_Line(" "); Put_Line("Boundary values computed from:"); Put_Line("ub(x,y,z) = x^4+2*y^4+3*z^4+4*x*y*z+5*x*y+z+2"); Put_Line("On cube -1 To 1 in three dimensions"); Put_Line(" "); Put_Line("xmin="&Real'Image(xmin)&", xmax="&Real'Image(xmax)); Put_Line("ymin="&Real'Image(ymin)&", ymax="&Real'Image(ymax)); Put_Line("zmin="&Real'Image(zmin)&", zmax="&Real'Image(zmax)); Make_Xnames; Start_time := Seconds(Clock); hx := (xmax-xmin)/real(nx-1); Put_Line("nx="&Integer'Image(nx)&" points for numeric derivative, step="& Real'Image(hx)); Put_Line("x grid"); 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)&" points for numeric derivative, step="& Real'Image(hy)); Put_Line("y grid"); for ii in 1..ny loop yg(ii) := ymin+Real(ii-1)*hy; end loop; hz := (zmax-zmin)/real(nz-1); Put_Line("nz="&Integer'Image(ny)&" points for numeric derivative, step="& Real'Image(hz)); Put_Line("z grid"); for iii in 1..nz loop zg(iii) := zmin+Real(iii-1)*hz; end loop; for i in 1..nx loop for ii in 1..ny loop for iii in 1..nz loop Ua(i,ii,iii) := ub(xg(i),yg(ii),zg(iii)); Put("xg="); Put(xg(i),3,5,0); Put(", yg="); Put(yg(ii),3,5,0); Put(", zg="); Put(zg(iii),3,5,0); Put(", Ua("&Integer'Image(i)&","&Integer'Image(ii)&","& ","&Integer'Image(iii)&")="); Put(Ua(i,ii,iii),3,5,0); Put(", f(x)="); Put(f(xg(i),yg(ii),zg(iii)),3,5,0); New_Line; if i=1 or i=nx or ii=1 or ii=ny or iii=1 or iii=nz then ug(i,ii,iii) := ub(xg(i),yg(ii),zg(iii)); else ug(i,ii,iii) := 0.0; end if; end loop; end loop; end loop; New_Line; Build_Vars; Compute_Discrete_Derivatives; Build_Equations; -- discrete Put_Line("unlin="&Integer'Image(Unlin)); Print_Equations(dof, unlin, Var1, Var2, Var3, Vari1, Vari2); Print_AY; Fit_Init(4,3); Fit_Bpoly(Nx, Ny, Nz, Xg, Yg, Zg, Ug); -- compute least sqaure fit for i in 2..nx-1 loop for ii in 2..ny-1 loop for iii in 2..nz-1 loop -- set up least square guess Xs(S(i,ii,iii)) := Eval_Poly(xg(i),yg(ii),zg(iii)); end loop; end loop; end loop; Put_Line("test, giving least square fit of boundary as initial guess"); Check_Equations(dof, unlin, A, Yy, Var1, Var2, Var3, Vari1, Vari2, Xs); Simeq_Newton5(dof, unlin, A, Yy, Var1, Var2, Var3, Vari1, Vari2, Xs, 20, 0.001, 1); Put_Line("observe convergence"); New_Line; for i in 2..nx-1 loop for ii in 2..ny-1 loop for iii in 2..nz-1 loop Xs(S(i,ii,iii)) := uA(i,ii,iii); -- set up exact solution end loop; end loop; end loop; Put_Line("test, giving exact solution"); Check_Equations(dof, unlin, A, Yy, Var1, Var2, Var3, Vari1, Vari2, Xs); Simeq_Newton5(dof, unlin, A, Yy, Var1, Var2, Var3, Vari1, Vari2, Xs, 20, 0.001, 1); New_Line; for i in 2..nx-1 loop for ii in 2..ny-1 loop for iii in 2..nz-1 loop Put("expected="); Put(Ua(i,ii,iii),8,3,0); Put(" got="); Put(Xs(S(i,ii,iii)),8,3,0); New_Line; end loop; end loop; end loop; New_Line; -- guess at initial conditions, boundary not included for I in 1..dof loop Xs(i) := 1.0; -- 5.0 did not converge end loop; Put_Line("initial guess, all "&Real'Image(Xs(1))); Simeq_Newton5(dof, unlin, A, Yy, Var1, Var2, Var3, Vari1, Vari2, Xs, 30, 0.0001, 1); for i in 2..nx-1 loop for ii in 2..ny-1 loop for iii in 2..nz-1 loop ug(i,ii,iii) := Xs(S(i,ii,iii)); -- ug is now boundary with -- computed solution end loop; end loop; end loop; avgerr := 0.0; maxerr := 0.0; Put_Line("U computed from nonlinear equations, Ua analytic, error "); for i in 1..nx loop for ii in 1..ny loop for iii in 1..nz loop err := abs(ug(i,ii,iii)-Ua(i,ii,iii)); if err>maxerr then maxerr := err; end if; avgerr := avgerr + err; Put("ug("&Integer'Image(i)&","&Integer'Image(ii)&","& Integer'Image(iii)&")="); Put(ug(i,ii,iii),4,8,0); Put(", Ua="); Put(Ua(i,ii,iii),4,8,0); Put(", err="); Put(ug(i,ii,iii)-Ua(i,ii,iii),4,8,0); New_Line; end loop; end loop; end loop; Put_Line(" maxerr="& Real'Image(maxerr)&", avgerr="& Real'image(avgerr/Real(dof))); New_Line; Put_line("check_soln on computed solution against PDE"); check_soln(ug); -- for use when analytic solution, Ua, is not known New_Line; Put_line("just checking check_soln when given correct solution"); check_soln(Ua); -- checking on checker New_Line; End_Time := Seconds(Clock); Put_Line("pde_nl33.adb finished in "&Duration'Image(End_Time-Start_Time)& " seconds"); end pde_nl33;