-- pde_nl22.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: -- -- a1(x,y)*uxx(x,y)*uyy(x,y) + b1(x,y)*u(x,y)*uxx(x,y) + -- c1(x,y)*u(x,y)uyy(x,y) = f(x,y) -- -- f(x,y) = 0.5*exp(x/2.0)*exp(y)*(6.0*x+6.0*y)* -- (12.0*y+8.0*x)+0.7/(x*x*y*y+0.5)* -- (x*x*x+2.0*y*y*y+3.0*x*x*y+4.0*x*y*y+ -- 5.0*x*y+6.0*x+7.0*y+8.0)*(6.0*x+6.0*y)+ -- (8.0-2.0*exp(x)-2.0*exp(0.5*y))* -- (x*x*x+2.0*y*y*y+3.0*x*x*y+4.0*x*y*y+ -- 5.0*x*y+6.0*x+7.0*y+8.0)*(12.0*y+8.0*x) -- -- a1(x,y) = exp(x/2)*exp(y)/2 -- b1(x,y) = 0.7/(x^2*y^2+0.5) -- c1(x,y) = 2*(4-exp(x)-exp(y/2)) -- -- Boundary conditions computed from -- ub(x,y) = x^4 + 2 y^4 + 3 x^3 y + 4 x y^3 + -- 5 x^2 y^2 + 6 x^2 + 7 y^2 + 8 -- -- Geometry of problem solution design, dof is linear degrees of freedom -- nx -- | 1 | 2 | 3 | 4 | 5 | -> X, cx(i,j), cxx(i,j) -- -----+-------+-------+-------+-------+-------+ i, j subscripts -- |(1,1) |(2,1) |(3,1) |(4,1) |(5,1) | -- 1 | | | | | | -- | B| B| B| B| B| -- -----+-------+-------+-------+-------+-------+ -- |(1,2) |(2,2) |(3,2) |(4,2) |(5,2) | <- ug subscript -- 2 | | X11 | X21 | X31 | | B are known -- | B| k=1| k=4| k=7| B| boundary values -- -----+-------+-------+-------+-------+-------+ -- |(1,3) |(2,3) |(3,3) |(4,3) |(5,3) | -- 3 | | X12 | X22 | X32 | | <- unknown linear -- | B| k=2| k=5| k=8| B| variables -- -----+-------+-------+-------+-------+-------+ -- |(1,4) |(2,4) |(3,4) |(4,4) |(5,4) | -- 4 | | X13 | X23 | X33 | | -- | B| k=3| k=6| k=9| B| <- k is variable subscript -- -----+-------+-------+-------+-------+-------+ -- ny |(1,5) |(2,5) |(3,5) |(4,5) |(5,5) | dof = (nx-2)*(ny-2) -- 5 | | | | | | nlin = dof*dof -- | B| B| B| B| B| ntot = dof+nlin -- -----+-------+-------+-------+-------+-------+ -- | -- V -- Y, cy, cyy k = S(i,ii), S(j,ii), S(i,jj), S(j,ii) -- ii,jj subscripts S(i,ii) = (i-2)*(ny-2)+(j-2)+1 -- this generalizes to 4 dimensions -- -- nonlinear terms k=10..18 X11*X11 .. X11*X33 nonlinear k = -- k=19..27 X12*X11 .. X12*X33 (k1-1)*dof+(k2-1)+1 -- k=28..36 X13*X11 .. X13*X33 for k1,k2 linear k -- k=35..45 X21*X11 .. X21*X33 dof=9 this example -- ... -- k=82..90 X33*X11 .. X33*X33 yes, redundancy for ease -- -- the nonlinear system of equations to solve is A X = Y -- dof equations in dof linear terms and nlin nonlinear terms -- -- | A(1,1) A(1,2) ... A(1,9) A(1,10) .. A(1,90) | |X11 | |Y(1)| -- | A(2,1) A(2,2) ... A(2,9) A(2,10) .. A(2,90) | |X12 | |Y(2)| -- | |*|... |=|... | -- | A(9,1) A(9,2) ... A(9,9) A(9,10) .. A(9,90) | |X33 | |Y(9)| -- |X11*X11| -- |... | -- |X33*X33| -- for discretization using 5 points, 4th order -- cx is the discretization coefficients, first subscript is where computed -- the derivatives at (2,2) (2,3) (3,2) are -- -- ux(2,2) = cx(2,1)*ug(1,2) + cx(2,2)*ug(2,2) + cx(2,3)*ug(3,2) + -- cx(2,4)*ug(4,2) + cx(2,5)*ug(5,2) note ug(1,2) and ug(5,2) known -- -- uy(2,2) = cy(2,1)*ug(2,1) + cy(2,2)*ug(2,2) + cy(2,3)*ug(2,3) + -- cy(2,4)*ug(2,4) + cy(2,5)*ug(2,5) note ug(2,1) and ug(2,5) known -- -- -- ux(3,2) = cx(3,1)*ug(1,2) + cx(3,2)*ug(2,2) + cx(3,3)*ug(3,2) + -- cx(3,4)*ug(4,2) + cx(3,5)*ug(5,2) note ug(1,2) and ug(1,5) known -- -- uy(3,2) = cy(2,1)*ug(3,1) + cy(2,2)*ug(3,2) + cy(2,3)*ug(3,3) + -- cy(2,4)*ug(3,4) + cy(2,5)*ug(3,5) note ug(3,1) and ug(3,5) known -- -- -- ux(2,3) = cx(2,1)*ug(1,3) + cx(2,2)*ug(2,3) + cx(2,3)*ug(3,3) + -- cx(2,4)*ug(4,3) + cx(2,5)*ug(5,3) note ug(1,3) and ug(5,3) known -- -- uy(2,3) = cy(3,1)*ug(2,1) + cy(3,2)*ug(2,2) + cy(3,3)*ug(2,3) + -- cy(3,4)*ug(2,4) + cy(3,5)*ug(2,5) note ug(2,1) and ug(2,5) known -- -- same subscripting for uxx using cxx, uyy using cyy -- -- OK, I know this is grossly elaborate, yet this is a teaching example. -- 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; procedure pde_nl22 is nx : Integer := 5; -- number of coefficients for discrete derivative ny : Integer := 5; dof : Integer := (nx-2)*(ny-2); -- number of equations and -- number of linear unknowns nlin : Integer := dof*dof; -- number of nonlinear terms, -- specific this PDE, nx-2 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 yg : Real_Vector(1..ny); -- y grid ug : Real_Matrix(1..nx,1..ny); -- computed solution at xy 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_matrix(1..nx,1..ny); -- known includes boundary cx : Real_Matrix(1..nx,1..nx); -- numerical first derivative cxx : Real_Matrix(1..nx,1..nx); -- numerical second derivative hx : real; cy : Real_Matrix(1..ny,1..ny); -- numerical first derivative cyy : Real_Matrix(1..ny,1..ny); -- numerical second derivative hy : real; err, avgerr, Maxerr, x , y: real; eps : real := 1.0e-4; xmin : Real := -1.0; xmax : Real := 1.0; ymin : Real := -1.0; ymax : Real := 1.0; Xnames : array(1..nx*ny) of String(1..3); -- nx<=9, ny<=9 -- e.g. ("X11", "X12", "X13", "X21", "X22", "X23", "X31", "X32", "X33"); Start_Time, End_Time : Duration; Debug : Boolean := False; function S(i:Integer; ii:Integer) return Integer is begin return (i-2)*(Ny-2)+(ii-2)+1; -- make single subscript end S; -- definition of PDE to be solved function a1(x:Real; y:Real) return real is begin return exp(x/2.0)*exp(y)/2.0; end a1; function b1(x:Real; y:Real) return real is begin return 0.7/(x*x*y*y+0.5); end b1; function c1(x:Real; y:Real) return real is begin return (4.0 - exp(x) - exp(y/2.0))*2.0; end c1; function f(x:Real; y:Real) return real is -- forcing function begin return 0.5*exp(x/2.0)*exp(y)*(6.0*x+6.0*y)* (12.0*y+8.0*x)+0.7/(x*x*y*y+0.5)* (x*x*x+2.0*y*y*y+3.0*x*x*y+4.0*x*y*y+ 5.0*x*y+6.0*x+7.0*y+8.0)*(6.0*x+6.0*y)+ (8.0-2.0*exp(x)-2.0*exp(0.5*y))* (x*x*x+2.0*y*y*y+3.0*x*x*y+4.0*x*y*y+ 5.0*x*y+6.0*x+7.0*y+8.0)*(12.0*y+8.0*x); end f; function ub(x:Real; y:Real) return real is -- analytic solution and boundary begin return x*x*x + 2.0*y*y*y + 3.0*x*x*y + 4.0*x*y*y + 5.0*x*y + 6.0*x + 7.0*y + 8.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 -- based on PDE k : Integer; 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; k := dof+1; for i in 1..dof loop for j in 1..dof loop -- set var1, var2 for nonlinear terms var1(k) := i; -- x1^2 x2^2 x3^2 var2(k) := j; var3(k) := 0; vari1(k) := 0; vari2(k) := 0; k := k+1; end loop; end loop; end Build_Vars; procedure Compute_Discrete_Derivatives is T : Real_Vector(1..nx+ny); begin -- first subscript is location i for cx, cxx ii for cy, cyy -- second subscript is point j for cx, cxx jj for cy, cyy 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(2, nx, i, hx, T); for j in 1..nx loop cxx(i,j) := T(j); if Debug then Put_Line("cxx("&Integer'Image(i)&","&Integer'Image(j)&")="& Real'Image(cxx(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(2, ny, ii, hy, T); for jj in 1..ny loop cyy(ii,jj) := T(jj); if Debug then Put_Line("cyy("&Integer'Image(ii)&","&Integer'Image(jj)&")="& Real'Image(cyy(ii,jj))); end if; -- debug end loop; end loop; end Compute_Discrete_Derivatives; procedure Build_Equations is ia, k : Integer; a1xy, b1xy, C1xy : Real; 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 -- a1(x,y)*uxx(x,y)*uyy(x,y) + b1(x,y)*u(x,y)*uxx(x,y) + -- c1(x,y)*u(x,y)uyy(x,y) = f(x,y) -- 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, X11..X33 linear terms y := yg(ii); ia := S(i,ii); -- equation index Yy(ia) := f(x,y); -- RHS a1xy := a1(x,y); b1xy := b1(x,y); c1xy := c1(x,y); -- term + a1(x,y)*uxx(x,y)*uyy(x,y) x at i, y at ii for j in 1..nx loop -- ug(i,j) cx(i,j) if j=1 or j=nx then -- boundary in x direction for jj in 1..ny loop -- ug(ii,jj) cy(ii,jj) if jj=1 or jj=ny then -- boundary in y direction and y direction if Debug then Put_Line("a1 i="&Integer'Image(i)&", j="&Integer'Image(j)& ", ii="&Integer'Image(ii)&", jj="&Integer'Image(jj)); Put_Line("cxx(i,j)="&Real'Image(cxx(i,j))& ", ug(i,j)="&Real'Image(ug(i,j))); Put_Line("cyy(ii,jj)="&Real'Image(cyy(ii,jj))& ", ug(ii,jj)="&Real'Image(ug(ii,jj))); end if; Yy(ia) := Yy(ia) - a1xy*(cxx(i,j)*ug(j,ii))* (cyy(ii,jj)*ug(i,jj)); else -- no boundary in y direction, boundary in x direction k := S(i,jj); -- linear u_i,jj if Debug then Put_Line("a1 i="&Integer'Image(i)&", j="&Integer'Image(j)& ", ii="&Integer'Image(ii)&", jj="&Integer'Image(jj)); Put_Line("S(i,jj)="&Integer'Image(S(i,jj))& ", k="&Integer'Image(k)); Put_Line("cxx(i,j)="&Real'Image(cxx(i,j))& ", ug(j,ii)="&Real'Image(ug(j,ii))); end if; -- debug A(ia,k) := A(ia,k) + a1xy*cyy(ii,jj)*(cxx(i,j)*ug(j,ii)); end if; end loop; else -- no boundary in x direction for jj in 1..ny loop -- ug(ii,jj) cy(ii,jj) if jj=1 or jj=ny then -- boundary in y direction k := S(j,ii); -- linear u_j,ii if Debug then Put_Line("a1 i="&Integer'Image(i)&", j="&Integer'Image(j)& ", ii="&Integer'Image(ii)&", jj="&Integer'Image(jj)); Put_Line("S(j,ii)="&Integer'Image(S(j,ii))& ", k="&Integer'Image(k)); Put_Line("cyy(ii,jj)="&Real'Image(cyy(ii,jj))& ", ug(i,jj)="&Real'Image(ug(i,jj))); end if; -- debug A(ia,k) := A(ia,k) + a1xy*cxx(i,j)*(cyy(ii,jj)*ug(i,jj)); else -- no boundary in y direction and x direction k := dof*S(j,ii) + S(i,jj); -- nonlinear u_j,ii*u_i,jj if Debug then Put_Line("a1 i="&Integer'Image(i)&", j="&Integer'Image(j)& ", ii="&Integer'Image(ii)&", jj="&Integer'Image(jj)); Put_Line("S(j,ii)="&Integer'Image(S(j,ii))& ", S(i,jj)="&Integer'Image(S(i,jj))& ", dof="&Integer'Image(Dof)& ", k="&Integer'Image(k)&" nl"); end if; -- debug A(ia,k) := A(ia,k) + a1xy*cxx(i,j)*cyy(ii,jj); end if; end loop; end if; end loop; -- nx for a1 -- term + b1(x,y)*u(x,y)*uxx(x,y) x at i, y at ii for j in 1..nx loop -- xg(i,j) cx(i,j) if j=1 or j=nx then -- boundary in x direction k := S(i,ii); -- linear at u if Debug then Put_Line("b1 i="&Integer'Image(i)&", j="&Integer'Image(j)& ", ii="&Integer'Image(ii)&", jj= none"); Put_Line("S(i,ii)="&Integer'Image(S(i,ii))& ", k="&Integer'Image(k)); Put_Line("cxx(i,j)="&Real'Image(cxx(i,j))& ", ug(i,j)="&Real'Image(ug(i,j))); end if; -- debug A(ia,k) := A(ia,k) + b1xy*(cxx(i,j)*ug(j,ii)); else -- no boundary in x direction k := dof*S(i,ii) + S(j,ii); -- nonlinear u*u_i,j if Debug then Put_Line("b1 i="&Integer'Image(i)&", j="&Integer'Image(j)& ", ii="&Integer'Image(ii)&", jj=none"); Put_Line("S(i,ii)="&Integer'Image(S(i,ii))& ", S(j,ii)="&Integer'Image(S(j,ii))& ", dof="&Integer'Image(Dof)& ", k="&Integer'Image(k)&" nl"); end if; -- debug A(ia,k) := A(ia,k) + b1xy*cxx(i,j); end if; end loop; -- nx for b1 -- term + c1(x,y)*u(x,y)uyy(x,y) x at i, y at ii for jj in 1..ny loop -- ug(ii,jj) cy(ii,jj) if jj=1 or jj=ny then -- boundary in y direction k := S(i,ii); if Debug then Put_Line("c1 i="&Integer'Image(i)&", j= none"& ", ii="&Integer'Image(ii)&", jj="&Integer'Image(jj)); Put_Line("S(i,ii)="&Integer'Image(S(i,ii))& ", k="&Integer'Image(k)); Put_Line("cyy(ii,jj)="&Real'Image(cyy(ii,jj))& ", ug(ii,jj)="&Real'Image(ug(ii,jj))); end if; -- debug A(ia,k) := A(ia,k) + c1xy*(cyy(ii,jj)*ug(i,jj)); else -- no boundary in y direction k := dof*S(i,ii) + S(i,jj); -- nonlinear u*u_i,jj if Debug then Put_Line("c1 i="&Integer'Image(i)&", j=none"& ", ii="&Integer'Image(ii)&", jj="&Integer'Image(jj)); Put_Line("S(i,ii)="&Integer'Image(S(i,ii))& ", S(i,jj)="&Integer'Image(S(i,jj))& ", dof="&Integer'Image(Dof)& ", k="&Integer'Image(k)&" nl"); end if; -- debug A(ia,k) := A(ia,k) + c1xy*cyy(ii,jj); end if; end loop; -- ny for c1 end loop; -- ii Y end loop; -- i X 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..ntot 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(Ux : Real_Matrix) is -- check when solution unknown smaxerr, err, x, y, uxx, uyy : Real; begin -- a1(x,y)*uxx(x,y)*uyy(x,y) + b1(x,y)*u(x,y)*uxx(x,y) + -- c1(x,y)*u(x,y)uyy(x,y) = f(x,y) 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); err := -f(x,y); -- add other side of equation terms uxx := 0.0; for j in 1..nx loop -- compute numerical derivative at (i,ii) uxx := uxx + cxx(i,j)*Ux(j,ii); end loop; uyy := 0.0; for jj in 1..ny loop uyy := uyy + cyy(ii,jj)*Ux(i,jj); end loop; err := err + a1(x,y)*uxx*uyy; err := err + b1(x,y)*Ux(i,ii)*uxx; err := err + c1(x,y)*Ux(i,ii)*uyy; if abs(err) > smaxerr then smaxerr := abs(err); end if; 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; nlin : 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+nlin 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+nlin then Put_Line("+"); else New_Line; end if; end loop; Put_Line(" = Y(i)"); New_Line; end Print_Equations; procedure Check_Equations(n : integer; nlin : 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+nlin 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+nlin 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 Xnames(S(i,ii)) := "X"&Integer'Image(i-1)(2..2)& Integer'Image(ii-1)(2..2); end loop; end loop; end Make_Xnames; begin Put_Line("pde_nl22.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("a1(x,y)*uxx(x,y)*uyy(x,y) + b1(x,y)*u(x,y)*uxx(x,y) +"); Put_Line("c1(x,y)*u(x,y)uyy(x,y) = f(x,y) "); Put_Line(" "); Put_Line("f(x,y)=0.5*exp(x/2.0)*exp(y)*(6.0*x+6.0*y)*"); Put_Line(" (12.0*y+8.0*x)+0.7/(x*x*y*y+0.5)*"); Put_Line(" (x*x*x+2.0*y*y*y+3.0*x*x*y+4.0*x*y*y+"); Put_Line(" 5.0*x*y+6.0*x+7.0*y+8.0)*(6.0*x+6.0*y)+"); Put_Line(" (8.0-2.0*exp(x)-2.0*exp(0.5*y))*"); Put_Line(" (x*x*x+2.0*y*y*y+3.0*x*x*y+4.0*x*y*y+"); Put_Line(" 5.0*x*y+6.0*x+7.0*y+8.0)*(12.0*y+8.0*x);"); Put_Line(" "); Put_Line("a1(x,y) = exp(x/2)*exp(y)/2"); Put_Line("b1(x,y) = 0.7/(x^2*y^2+0.5)"); Put_Line("c1(x,y) = 2*(4-exp(x)-exp(y/2))"); Put_Line(" "); Put_Line("Boundary values computed from:"); Put_Line("ub(x,y)= x^4 + 2 y^4 + 3 x^3 y + 4 x y^3 +"); Put_Line(" 5 x^2 y^2 + 6 x^2 + 7 y^2 + 8"); Put_Line(" "); Put_Line("xmin="&Real'Image(xmin)&", xmax="&Real'Image(xmax)); Put_Line("ymin="&Real'Image(ymin)&", ymax="&Real'Image(ymax)); 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; for i in 1..nx loop for ii in 1..ny loop Ua(i,ii) := ub(xg(i),yg(ii)); Put("xg="); Put(xg(i),3,5,0); Put(", yg="); Put(yg(ii),3,5,0); Put(", Ua("&Integer'Image(i)&","&Integer'Image(ii)&")="); Put(Ua(i,ii),3,5,0); Put(", f(x)="); Put(f(xg(i),yg(ii)),3,5,0); New_Line; if i=1 or i=nx or ii=1 or ii=ny then ug(i,ii) := ub(xg(i),yg(ii)); else ug(i,ii) := 0.0; end if; end loop; end loop; New_Line; Build_Vars; Compute_Discrete_Derivatives; Build_Equations; -- discrete Print_AY; Print_Equations(dof, nlin, Var1, Var2, Var3, Vari1, Vari2); for i in 2..nx-1 loop for ii in 2..ny-1 loop Xs(S(i,ii)) := uA(i,ii); -- set up exact solution end loop; end loop; Put_Line("test, giving exact solution"); Check_Equations(dof, nlin, A, Yy, Var1, Var2, Var3, Vari1, Vari2, Xs); Simeq_Newton5(dof, nlin, 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 Put("expected="); Put(Ua(i,ii),8,3,0); Put(" got="); Put(Xs(S(i,ii)),8,3,0); New_Line; end loop; end loop; New_Line; -- guess at initial conditions, boundary not included for I in 1..dof loop Xs(i) := 5.0; -- 1.0 did not converge end loop; Put_Line("initial guess, all "&Real'Image(Xs(1))); Simeq_Newton5(dof, nlin, A, Yy, Var1, Var2, Var3, Vari1, Vari2, Xs, 30, 0.001, 1); for i in 2..nx-1 loop for ii in 2..ny-1 loop ug(i,ii) := Xs(S(i,ii)); -- ug is now boundary with computed solution 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 err := abs(ug(i,ii)-Ua(i,ii)); if err>maxerr then maxerr := err; end if; avgerr := avgerr + err; Put("ug("&Integer'Image(i)&","&Integer'Image(ii)&")="); Put(ug(i,ii),4,8,0); Put(", Ua="); Put(Ua(i,ii),4,8,0); Put(", err="); Put(ug(i,ii)-Ua(i,ii),4,8,0); New_Line; 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_nl22.adb finished in "&Duration'Image(End_Time-Start_Time)& " seconds"); end pde_nl22;