program ConventionalPoisson4thOrder_Example_2;
{$APPTYPE CONSOLE}
uses SysUtils, DateUtils;
type matrix        = array of array of Extended;
     rightfunction = function (const x, y : Extended) : Extended;
     boundary      = function (const x : Extended) : Extended;

procedure fullPoisson (const n, m                   : Integer;
                       const alpha, beta, eps       : Extended;
                       const f                      : rightfunction;
                       const phi1, phi2, phi3, phi4 : boundary;
                       out u                        : matrix;
                       out st                       : Integer);
{------------------------------------------------------------------------------}
{                                                                              }
{  The procedure fullPoisson solves the second order partial differential      }
{  equation u|xx + u|yy = f(x,y), where u=u(x,y), u|xx denotes the second      }
{  order partial derivative of u with respect to x, 0 <= x <= alpha,           }
{  0 <= y <= beta, with the boundary conditions u = phi1(y) for x = 0,         }
{  u = phi2(x) for y = 0, u = phi3(y) for x = alpha, u = phi4(x) for y = beta, }
{  where phi1(0) = phi2(0), phi2(alpha) = phi3(0), phi3(beta) = phi4(alpha),   }
{  phi4(0) = phi1(beta). The problem is solved by a difference method of 4th   }
{  order.                                                                      }
{                                                                              }
{  Data:                                                                       }
{    n           - the number of subintervals in the interval [0,alpha],       }
{    m           - the number of subintervals in the interval [0,beta],        }
{    alpha       - a maximum value of x,                                       }
{    beta        - a maximum value of y,                                       }
{    eps         - the maximum difference for which the values of phi1(0) and  }
{                  phi2(0), phi2(alpha) and phi3(0), phi3(beta) and            }
{                  phi4(alpha), phi4(0) and phi1(beta) can be considered as    }
{                  equal,                                                      }
{    f           - a Delphi Pascal function which calculates the value of      }
{                  function f = f(x,y) occurring in the differential equation, }
{    phi1, phi2,                                                               }
{    phi3, phi4  - Delphi Pascal functions which calculate the values of       }
{                  functions phi1 = phi1(y), phi2 = phi2(x), phi3 = phi3(y)    }
{                  and phi4 = phi4(x) occurring in the boundary conditions.    }
{  Result:                                                                     }
{    u - an array containing the solution (the element u[i,j] contains the     }
{        solution at the point (ih, jk), where i = 0, 1, ... , n,              }
{        j = 0, 1, ... , m, h = alpha/n, k = beta/m.                           }
{  Other parameters:                                                           }
{    st   - a variable which within the procedure diffmethod is assigned the   }
{           value of:                                                          }
{             1, if n < 2 or m < 2,                                            }
{             2, if alpha <= 0 or beta <= 0,                         }
{             4, if at least one of the conditions given above in the          }
{                description of eps is not fulfilled (that means the boundary  }
{                conditions are wrong),                                        }
{             5, if the matrix of adequate linear system is singular,          }
{             0, otherwise.                                                    }
{           Note: If st<>0, then the elements of array u are not calculated.   }
{  Unlocal identifiers:                                                        }
{    matrix        - a type identifier of interval dynamical array defined as  }
{                    follows:                                                  }
{                       type matrix = array of array of Extended;              }
{                    which corresponds to the array [0..n, 0..m],              }
{    rightfunction - a procedural-type identifier defined as follows           }
{                       type rightfunction = function (const x, y : Extended)  }
{                                                        : Extended;           }
{    boundary      - a procedural-type identifier defined as follows           }
{                       type boundary = function (const X : Extended)          }
{                                                   : Extended;                }
{  Note:                                                                       }
{    Within the procedure it is assigned a memory for the u parameter          }
{    (dynamical array). After using the procedure you have to release this     }
{    memory by calling the procedure Finalize(x), where x denotes your         }
{    argument corresponding to the u parameter.                                }
{                                                                              }
{------------------------------------------------------------------------------}
var i, j, jh, j1, k, kh, l, lh, l1, l2, n1, n2, p, q, rh   : Integer;
    af, bf, cf, h1, h2, h4, ihh, iminus1hh, iplus1hh, jkk,
    jminus1kk, jplus1kk, k1, k2, k4, max, pq, s, s1, s2, z : Extended;
    a1, b1, x                                              : array of Extended;
    r                                                      : array of Integer;
function boundconds (const b1, b2, eps : Extended) : Integer;
begin
  Result:=0;
  if (b1<>0) and (b2<>0)
    then if Abs(b1-b2)/Abs(b1)>=eps
           then Result:=4
           else
    else if b1=0
            then if Abs(b2)>=eps
                   then Result:=4
                   else
            else if b2=0
                   then if Abs(b1)>=eps
                          then Result:=4;
end {boundconds};
begin
  st:=0;
  if (n<2) or (m<2)
    then st:=1
    else if (alpha<=0) or (beta<=0)
           then st:=2;
  if st=0
    then
      begin
        s:=phi1(0);
        s1:=phi2(0);
        st:=boundconds(s,s1,eps);
        if st=0
          then begin
                 s:=phi2(alpha);
                 s1:=phi3(0);
                 st:=boundconds(s,s1,eps);
                 if st=0
                   then begin
                          s:=phi3(beta);
                          s1:=phi4(alpha);
                          st:=boundconds(s,s1,eps);
                          if st=0
                            then begin
                                   s:=phi4(0);
                                   s1:=phi1(beta);
                                   st:=boundconds(s,s1,eps)
                                 end
                        end
               end
      end;
  if st=0
    then begin
           h1:=alpha/n;
           k1:=beta/m;
           h2:=h1*h1;
           k2:=k1*k1;
           h4:=h2*h2;
           k4:=k2*k2;
           af:=h2+k2;
           bf:=2*(5*h2-k2);
           cf:=2*(5*k2-h2);
           n1:=(n-1)*(m-1);
           n2:=n1+1;
           p:=n2;
           SetLength (r,n1+1);
           for i:=1 to n2 do
             r[i-1]:=0;
           k:=0;
           j:=0;
           SetLength (a1,n1+1);
           SetLength (b1,n1+1);
           SetLength (x,(n*m-n-m+4)*(n*m-n-m+4) div 4);
           repeat
             k:=k+1;
             for i:=1 to n1 do
               a1[i-1]:=0;
             j:=j+1;
             i:=Trunc((k-1)/(m-1))+1;
             l1:=(i-2)*(m-1)+j;
             l2:=l1+m-1;
             ihh:=i*h1;
             jkk:=j*k1;
             iminus1hh:=(i-1)*h1;
             iplus1hh:=(i+1)*h1;
             jminus1kk:=(j-1)*k1;
             jplus1kk:=(j+1)*k1;
             if i>1
               then begin
                      a1[l1-1]:=cf;
                      if j>1
                        then a1[l1-2]:=af;
                      if j<m-1
                        then a1[l1]:=af
                    end;
             a1[l2-1]:=-20*af;
             if j>1
               then a1[l2-2]:=bf;
             if j<m-1
               then a1[l2]:=bf;
             l1:=l2+m-1;
             if i<n-1
               then begin
                      a1[l1-1]:=cf;
                      if j>1
                        then a1[l1-2]:=af;
                      if j<m-1
                        then a1[l1]:=af
                    end;
             s:=f(iminus1hh,jkk);
             s1:=f(iplus1hh,jkk);
             s:=s+s1;
             s1:=f(ihh,jkk);
             s:=s+8*s1;
             s1:=f(ihh,jminus1kk);
             s:=s+s1;
             s1:=f(ihh,jplus1kk);
             s:=s+s1;
             s:=h2*k2*s;
             if i=1
               then begin
                      s1:=phi1(jminus1kk);
                      s2:=phi1(jplus1kk);
                      s1:=af*(s1+s2);
                      s2:=phi1(jkk);
                      s1:=s1+cf*s2;
                      s:=s-s1;
                      if j=1
                        then begin
                               s1:=phi2(iplus1hh);
                               s:=s-af*s1;
                               s1:=phi2(ihh);
                               s:=s-bf*s1
                             end;
                      if j=m-1
                        then begin
                               s1:=phi4(iplus1hh);
                               s:=s-af*s1;
                               s1:=phi4(ihh);
                               s:=s-bf*s1
                             end
                    end
               else if i=n-1
                      then begin
                             s1:=phi3(jminus1kk);
                             if st=0
                               then begin
                                      s2:=phi3(jplus1kk);
                                      s1:=af*(s1+s2);
                                      s2:=phi3(jkk);
                                      s1:=s1+cf*s2
                                    end;
                             s:=s-s1;
                             if j=1
                               then begin
                                      s1:=phi2(iminus1HH);
                                      s:=s-af*s1;
                                      s1:=phi2(ihh);
                                      s:=s-bf*s1
                                    end;
                             if j=m-1
                               then begin
                                      s1:=phi4(iminus1hh);
                                      s:=s-af*s1;
                                      s1:=phi4(ihh);
                                      s:=s-bf*s1
                                    end
                           end
                      else begin
                             if j=1
                               then begin
                                      s1:=phi2(iminus1hh);
                                      s2:=phi2(iplus1hh);
                                      s1:=s1+s2;
                                      s:=s-af*s1;
                                      s1:=phi2(ihh);
                                      s:=s-bf*s1
                                    end;
                             if j=m-1
                               then begin
                                      s1:=phi4(iminus1hh);
                                      s2:=phi4(iplus1hh);
                                      s1:=s1+s2;
                                      s:=s-af*s1;
                                      s1:=phi4(ihh);
                                      s:=s-bf*s1
                                    end
                           end;
             a1[n2-1]:=s;
             for i:=1 to n1 do
               begin
                 rh:=r[i-1];
                 if rh<>0
                   then b1[rh-1]:=a1[i-1]
               end;
             kh:=k-1;
             l:=0;
             max:=0;
             for j1:=1 to n2 do
               if r[j1-1]=0
                 then begin
                        s:=a1[j1-1];
                        l:=l+1;
                        q:=l;
                        for i:=1 to kh do
                          begin
                            s:=s-b1[i-1]*x[q-1];
                            q:=q+p
                          end;
                        a1[l-1]:=s;
                        if s<0
                          then s:=Abs(s);
                        if (j1<n2) and (s>max)
                          then begin
                                 max:=s;
                                 jh:=j1;
                                 lh:=l
                               end
                      end;
             if max=0
               then st:=5
               else begin
                      max:=1/a1[lh-1];
                      r[jh-1]:=k;
                      for i:=1 to p do
                        a1[i-1]:=max*a1[i-1];
                      jh:=0;
                      q:=0;
                      for j1:=1 to kh do
                        begin
                          s:=x[q+lh-1];
                          for i:=1 to p do
                            if i<>lh
                              then begin
                                     jh:=jh+1;
                                     x[jh-1]:=x[q+i-1]-s*a1[i-1]
                                   end;
                          q:=q+p
                        end;
                      for i:=1 to p do
                        if i<>lh
                          then begin
                                 jh:=jh+1;
                                 x[jh-1]:=a1[i-1]
                               end;
                      p:=p-1
                    end;
             if j=m-1
               then j:=0
           until (k=n1) or (st<>0);
           Finalize (a1);
           Finalize (b1);
           if st=0
             then begin
                    SetLength (u,n+1);
                    for i:=0 to n do
                      SetLength (u[i],m+1);
                    for k:=1 to n1 do
                      begin
                        rh:=r[k-1];
                        if rh<>k
                          then begin
                                 s:=x[k-1];
                                 x[k-1]:=x[rh-1];
                                 i:=r[rh-1];
                                 while i<>k do
                                   begin
                                     x[rh-1]:=x[i-1];
                                     r[rh-1]:=rh;
                                     rh:=i;
                                     i:=r[rh-1]
                                   end;
                                 x[rh-1]:=s;
                                 r[rh-1]:=rh
                               end
                      end;
                    for i:=1 to n-1 do
                      for j:=1 to m-1 do
                        u[i,j]:=x[(i-1)*(m-1)+j-1];
                    i:=0;
                    repeat
                      i:=i+1;
                      ihh:=i*h1;
                      u[i,0]:=phi2(ihh);
                      u[i,m]:=phi4(ihh)
                    until (i=n-1) or (st<>0);
                    j:=-1;
                    repeat
                      j:=j+1;
                      jkk:=j+k1;
                      u[0,j]:=phi1(jkk);
                      u[n,j]:=phi3(jkk)
                    until (j=m) or (st<>0)
                  end;
           Finalize (X);
           Finalize (r)
         end
end {fullPoisson};

//  The Delphi Pascal interval functions for the equation
//  u|xx + u|yy = -2*Pi^2*sin(Pi*x)*sin(Pi*y) with u = 0 on the boundary
//  (0<=x<=1, 0<=y<=1).

function f (const x, y : Extended) : Extended;
begin
  Result:=-2*Pi*Pi*Sin(Pi*x)*Sin(Pi*y)
end {f};

function phi1 (const y : Extended) : Extended;
begin
  Result:=0
end {phi1};

function phi2 (const x : Extended) : Extended;
begin
  Result:=0
end {phi2};

function phi3 (const y : Extended) : Extended;
begin
  Result:=0
end {phi3};

function phi4 (const x : Extended) : Extended;
begin
  Result:=0
end {phi4};

function modify_time (const time1, time2 : TDateTime) : string;
var code, days, hour_numbers : Integer;
    hours, time              : string;
begin
  time:=TimeToStr(time1-time2);
  days:=DaysBetween(time1,time2);
  if days>0
    then begin
           hours:=Copy(time,1,2);
           Val (hours, hour_numbers, code);
           hour_numbers:=24*days+hour_numbers;
           Str (hour_numbers, hours);
           Result:=hours+Copy(time,3,6)
         end
    else Result:=time
end {modify_time};

var i, imod, j, jmod, l, n, m, st, ui, uj : Integer;
    alpha, beta, eps, exact, h, k         : Extended;
    file_name, time                       : string;
    z, z1                                 : Char;
    OK, OK1                               : Boolean;
    plik                                  : TextFile;
    time1, time2                          : TDateTime;
    uu                                    : matrix;

begin
  Writeln ('SOLVING THE POISSON EQUATION WITH GIVEN BOUNDARY CONDITIONS');
  Writeln ('BY A DIFFERENCE METHOD OF FOURTH ORDER');
  alpha:=1;
  beta:=1;
  Writeln;
  Writeln ('Insert data');
  repeat
    OK:=false;
    try
      Write ('n = m = ');
      Readln (n);
      if (n>=20) and (n mod 10 = 0)
        then OK:=True
        else Writeln ('Error: the value of n = m must be an integer number '
                      +'equal to 20, 30, 40, ... .');
      except
        Writeln ('Error: the value of n = m must be an integer number equal to '
                 +'20, 30, 40, ... .')
    end
  until OK;
  m:=n;

// Don't change the value of eps!
  eps:=1.0E-16;

  Writeln;
  repeat
    Write ('Output of results (f - file, s - screen): ');
    Readln (z)
  until (z='f') or (z='F') or (z='s') or (z='S');
  if (z='f') or (z='F')
    then begin
           Writeln;
           repeat
             OK:=false;
             try
               Writeln ('Insert the name of file');
               Write ('(the extension .txt will be added automatically): ');
               Readln (file_name);
               file_name:=file_name+'.txt';
               OK1:=False;
               if FileExists(file_name)
                 then begin
                        Writeln ('The file '+file_name+' already exists. ');
                        repeat
                          Writeln ('Do you want to delete this file and create '
                                   +'a new one in its place?');
                          Write ('(y - yes, n - no): ');
                          Readln (z1)
                        until (z1='y') or (z1='Y') or (z1='n') or (z1='N');
                        if (z1='y') or (z1='Y')
                          then OK1:=True
                      end
                 else OK1:=True;
               if OK1
                 then begin
                        AssignFile (plik, file_name);
                        Rewrite (plik);
                        OK:=True;
                        Writeln (plik, 'SOLVING THE POISSON EQUATION WITH GIVEN'
                                 +' BOUNDARY CONDITIONS');
                        Writeln (plik, 'BY A DIFFERENCE METHOD OF FOURTH ORDER')
                      end
             except
               Writeln ('Error: incorrect name of the .txt file.')
             end
           until OK
         end;
  Writeln;
  Writeln ('Calculating ... Please, wait!');
  Writeln;
  time1:=Now;
  fullPoisson (n, m, alpha, beta, eps, f, phi1, phi2, phi3, phi4, uu, st);
  time2:=Now;
  time:=modify_time(time1,time2);
  Writeln ('status = ', st, ', time = ', time);
  if (z='f') or (z='F')
    then Writeln (plik, 'status = ', st, ', time = ', time);
  if (z='s') or (z='S')
    then begin
           Writeln ('  u - a solution obtained by interval method,');
           Writeln (' eu - the exact (approximate) solution');
           Readln
         end
    else begin
           Writeln (plik, '  u - a solution obtained by interval method,');
           Writeln (plik, ' eu - the exact (approximate) solution')
         end;
  if st=0
    then begin
           h:=alpha/n;
           k:=beta/m;
           l:=0;
           imod:=n div 10;
           jmod:=m div 10;
           ui:=n div 2;
           uj:=m div 2;
           Writeln;
           for j:=0 to m do
             if j mod jmod = 0
               then begin
                      exact:=Sin(Pi*ui*h)*Sin(Pi*j*k);
                      if (z='s') or (z='S')
                        then begin
                               Writeln (' u(', ui*h:3:1, ',', j*k:3:1, ') = ',
                                        uu[ui,j]:26);
                               Write ('eu(', ui*h:3:1, ',', j*k:3:1, ') = ',
                                      exact:26);
                               Writeln ('   difference =',
                                        Abs(uu[ui,j]-exact):11);
                               l:=l+1;
                               if l=5
                                 then begin
                                        l:=0;
                                        Readln
                                      end
                                 else Writeln
                             end
                        else begin
                               Writeln (plik, ' ');
                               Writeln (plik, ' u(', ui*h:3:1, ',', j*k:3:1,
                                        ') = ', uu[ui,j]:26);
                               Write (plik,'eu(', ui*h:3:1, ',', j*k:3:1,
                                      ') =  ', exact:26);
                               Writeln (plik, '   difference =',
                                        Abs(uu[ui,j]-exact):11);
                             end
                    end;
           for i:=0 to n do
             if i mod imod = 0
               then begin
                      exact:=Sin(Pi*i*h)*Sin(Pi*uj*k);
                      if (z='s') or (z='S')
                        then begin
                               Writeln (' u(', i*h:3:1, ',', uj*k:3:1, ') = ',
                                        uu[i,uj]:26);
                               Write ('eu(', i*h:3:1, ',', uj*k:3:1, ') = ',
                                      exact:26);
                               Writeln ('   difference =',
                                        Abs(uu[i,uj]-exact):11);
                               l:=l+1;
                               if l=5
                                 then begin
                                        l:=0;
                                        Readln
                                      end
                                 else Writeln
                             end
                        else begin
                               Writeln (plik, ' ');
                               Writeln (plik, ' u(', i*h:3:1, ',', uj*k:3:1,
                                        ') = ', uu[i,uj]:26);
                               Write (plik, 'eu(', i*h:3:1, ',', uj*k:3:1,
                                      ') =  ', exact:26);
                               Writeln (plik, '   difference =',
                                        Abs(uu[i,uj]-exact):11);
                             end
                    end;
           Finalize (uu)
         end;
  if (z='f') or (z='F')
    then CloseFile(plik);
  Writeln ('THE END');
  Readln
end.