(*==========================================================================*)
(*                                                                          *)
(*  PASCAL-XSC - MODUL  IVP_DATA                                    000221  *)
(*                                                                          *)
(*  MODUL DER KOEFFIZIENTENFUNKTIONEN DES BEHANDELTEN ANFANGSWERTPROBLEMS   *)
(*                                                                          *)
(*==========================================================================*)

MODULE ivp_data;

USE i_ari , c_ari , ci_ari , icf_ari , mp_ari , mpi_ari ,
    service ,  mps_aril , itaylor , mps_tayl;


(*--------------------------------------------------------------------------*)

GLOBAL PROCEDURE TAYKOEFF_P_I( BSPNR: INTEGER; ALPHA, BETA, GAMMA: REAL;
                               VAR TIB: MPMATRIX; X_0: REAL; PREC: INTEGER );

(*  IN DIESER PROZEDUR WERDEN DIE TAYLORKOEFFIZIENTEN          *)
(*  DER KOEFFIZIENTENFUNKTIONEN DER DGL BERECHNET              *)
(*  MUSS FUER JEDE DGL NEU PROGRAMMIERT UND COMPILIERT WERDEN  *)
(*  ALPHA, BETA UND GAMMA SIND FUNKTIONSPARAMETER              *)

VAR I , J                : INTEGER;
    MPI_PI , MPIA , MPIB : MPINTERVAL;
    MPA , MPB            : MULPREC;

BEGIN

   (*  BEISPIEL 1: Y" = ALPHA e^x Y + e^(-x) - ALPHA  *)

   IF BSPNR = 1 THEN BEGIN

      FOR I := 0 TO UB( TIB[0,*] ) DO BEGIN
          TIB[1,I]  := 0;
          TIB[0,I].PREC := PREC;
          TIB[0,I]  := EXP( TIB[0] , 1.0 , X_0 , I );
          TIB[-1,I].PREC := PREC;
          TIB[-1,I] := EXP( TIB[-1] , -1.0 , X_0 , I );
      END;
      FOR I := 0 TO UB( TIB[0,*] ) DO
          TIB[0,I] := ALPHA *< TIB[0,I];
      TIB[-1,0] := TIB[-1,0] - MPS( ALPHA );

   END;

   (*  BEISPIEL 2: Y" = - ( ALPHA + BETA * COS( GAMMA * X ) ) * Y  *)

   IF BSPNR = 2 THEN BEGIN

      FOR I := 0 TO UB( TIB[0,*] ) DO BEGIN
          TIB[1,I]  := 0;
          TIB[0,I].PREC := PREC;
          TIB[0,I] := COS( TIB[0] , GAMMA , X_0 , I );
      END;

      FOR I := 0 TO UB( TIB[0,*] ) DO
          TIB[0,I] := - BETA *< TIB[0,I];

      TIB[0,0] := TIB[0,0] - MPS( ALPHA );

   END;

   (*  BEISPIEL 3: Y" = 400 Y + 400 COS^2( PI*X ) + 2 * PI^2 COS( 2*PI*X )  *)
   (*                 = 400 Y + ( 200 + 2 * PI^2 ) COS( 2*PI*X ) + 200      *)

   IF BSPNR = 3 THEN BEGIN

      MPINIT( MPIA );
      MPINIT( MPIB );
      MPINIT( MPI_PI );
      MPIA := 1;
      MPI_PI := 4 * ARCTAN( MPIA );
      MPIB := SQR( 2 * MPI_PI );
      MPA := CONVERT( MPIB , PREC );
      MPIB := 200 + 2 * SQR( MPI_PI );
      MPB := CONVERT( MPIB , PREC );

      TIB[1,0]  := 0;
      TIB[0,0]  := 400;
      MPIB := COS( ( 2 * X_0 ) * MPI_PI );
      TIB[-1,0] := CONVERT( MPIB , PREC );
      TIB[1,1]  := 0;
      TIB[0,1]  := 0;
      MPIB := - ( 2 * MPI_PI ) * SIN( ( 2 * X_0 ) * MPI_PI );
      TIB[-1,1] := CONVERT( MPIB , PREC );
      FOR I := 2 TO UB( TIB[0,*] ) DO BEGIN
          TIB[1,I]  := 0;
          TIB[0,I]  := 0;
          TIB[-1,I] := ( - MPA *< TIB[-1,I-2] ) / ( I * ( I - 1 ) );
      END;
      FOR I := 0 TO UB( TIB[0,*] ) DO
          TIB[-1,I] := MPB *< TIB[-1,I];
      TIB[-1,0] := TIB[-1,0] + MPS( 200 );

   (* FOR I := 0 TO UB( TIB[0,*] ) DO MPSWRITE( 'TIB[-1,I]' , TIB[-1,I] ); *)

      MPFREE( MPIA );
      MPFREE( MPIB );
      MPFREE( MPI_PI );

   END;

   (*  BEISPIEL 4: Y^(4) = - 6 ALPHA / ( 1 + ALPHA + X ) Y^(3)    *)
   (*                      - 6 ALPHA^2 / ( 1 + ALPHA + X )^2 Y''  *)
   (*                      - BETA / ( 1 + ALPHA + X )^2 Y         *)

   IF BSPNR = 4 THEN BEGIN

      TIB[3,0] := - 1;
      TIB[2,0] := - 1;
      TIB[1,0] := 0;
      FOR I := 1 TO UB( TIB[0,*] ) DO BEGIN
          TIB[3,I]  := ROUND( - TIB[3,I-1] *> ALPHA , PREC );
          TIB[2,I]  := ROUND( TIB[3,I] *> ( I + 1 ) , PREC );
          TIB[1,I]  := 0;
      END;
      FOR I := 0 TO UB( TIB[0,*] ) DO BEGIN
          TIB[0,I] :=  ROUND( TIB[2,I] *> BETA , PREC );
          TIB[2,I] :=  ROUND( ( ( TIB[2,I] *> ALPHA ) *> ALPHA ) *> 6 , PREC );
          TIB[3,I] :=  ROUND( ( TIB[3,I] *> ALPHA ) *> 6 , PREC );
      END;

   END;

END;

(*--------------------------------------------------------------------------*)

GLOBAL PROCEDURE F_A( BSPNR: INTEGER; ALPHA, BETA, GAMMA, MITTE, RADIUS: REAL;
                      PHI: INTERVAL; POLYORD: INTEGER; TKOEFF: MPVECTOR;
                      VAR CIMIDFUNKT , CIFUNKT: CINTERVAL; VAR MINDIAM: REAL;
                      VAR ERR: INTEGER );

(*  FUNKTIONSWERTE MIT SUBTRAKTION EINES APPROXIMIERENDEN TAYLORPOLYNOMES  *)

(*  DEFINITION DER IM MODUL CAUCHYCO ZU UNTERSUCHENDEN FUNKTIONEN  *)
(*  MUSS FUER JEDE DGL NEU PROGRAMMIERT UND COMPILIERT WERDEN      *)
(*  (EVENTUELL MIT CHECK DES DEFINITIONSBEREICHS)                  *)

(*  EINGABEVARIABLEN:                                              *)
(*  ALPHA, BETA...: FUNKTIONSPARAMETER                             *)
(*  MITTE, RADIUS : DES KREISES, AUF DEM F AUSZUWERTEN IST         *)
(*  PHI           : WINKEL DES AKTUELLEN KREISSEGMENTS             *)
(*  POLYORD       : ORDNUNG DES APPROXIMIERENDEN TAYLORPOLYNOMS    *)
(*  TKOEFF        : TAYLORKOEFFIZIENTEN                            *)
(*  MINDIAM       : HALBER DURCHMESSER DER FUNKTIONSAUSWERTUNG     *)
(*                  OHNE SUBTRAKTION DES TAYLORPOLYNOMS            *)
(*                  (UNTERSCHRANKE DER APPROXIMATIONSGUETE)        *)

VAR I , J , L                      : INTEGER;
    PHI_MID , ONE , IPI            : INTERVAL;
    CIEIN , CIMID , PINT , PINTMID : CINTERVAL;
    VIC , VID                      : IVECTOR[0..POLYORD+2];

BEGIN

   CIEIN.RE := MITTE + RADIUS * COS( PHI );
   CIEIN.IM := RADIUS * SIN( PHI );
   PHI_MID := MID( PHI );
   CIMID.RE := MITTE + RADIUS * COS( PHI_MID );
   CIMID.IM := RADIUS * SIN( PHI_MID );

   (*  ZUERST FUNKTIONSWERT VON  F  BERECHNEN  *)

   (*  BEISPIEL 1: Y" = ALPHA e^x Y + e^(-x) - ALPHA  *)

   IF BSPNR = 1 THEN BEGIN
      (*  F( Z ) = ALPHA * EXP( Z )  *)
      CIFUNKT := ALPHA * EXP( CIEIN );
      CIMIDFUNKT := ALPHA * EXP( CIMID );
   END;

   (*  BEISPIEL 2: Y" = - ( ALPHA + BETA * COS( GAMMA * X ) ) * Y  *)

   IF BSPNR = 2 THEN BEGIN
      (*  F( Z ) = - ALPHA - BETA * COS( GAMMA*Z )  *)
      CIFUNKT := - ALPHA - BETA * COS( GAMMA * CIEIN );
      CIMIDFUNKT := - ALPHA - BETA * COS( GAMMA * CIMID );
   END;

   (*  BEISPIEL 3: Y" = 400 Y + 400 COS^2( PI*X ) + 2 * PI^2 COS( 2*PI*X )  *)
   (*                 = 400 Y + ( 200 + 2 * PI^2 ) COS( 2*PI*X ) + 200      *)

   IF BSPNR = 3 THEN BEGIN
      (*  F( Z ) = 200 + ( 200 + 2 * SQR( PI ) ) COS( 2*PI*X )  *)
      ONE := 1;
      IPI := 4 * ARCTAN( 1 );
      CIFUNKT := 200 + ( 200 + 2 * SQR( IPI ) ) * COS( 2 * IPI * CIEIN );
      CIMIDFUNKT := 200 + ( 200 + 2 * SQR( IPI ) ) * COS( 2 * IPI * CIMID );
   END;

   (*  BEISPIEL 4: Y^(4) = - 6 ALPHA / ( 1 + ALPHA * X ) Y^(3)    *)
   (*                      - 6 ALPHA^2 / ( 1 + ALPHA * X )^2 Y''  *)
   (*                      - BETA / ( 1 + ALPHA * X )^2 Y         *)

   IF BSPNR = 4 THEN BEGIN
      (*  F( Z ) = 1 / ( 1 + ALPHA * Z )  *)
      CIFUNKT := 1 / ( 1 + ALPHA * CIEIN );
      CIMIDFUNKT := 1 / ( 1 + ALPHA * CIMID );
   END;

   IF BSPNR = 41 THEN BEGIN
      (*  F( Z ) = 1 / ( 1 + ALPHA * Z )^2  *)
      CIFUNKT := SQR( 1 / ( 1 + ALPHA * CIEIN ) );
      CIMIDFUNKT := SQR( 1 / ( 1 + ALPHA * CIMID ) );
   END;

   (*  BEISPIEL ? MIT LOGARITHMUSFUNKTION  *)

   (*  F( Z ) = LN( Z )                       *)
   (*  DEFINITIONSBEREICH IST ZU ÜBERPRÜFEN!  *)
(* CIFUNKT := LNP( CIEIN , CIFUNKT , ERR );*)
(* CIMIDFUNKT := LNP( CIMID , CIMIDFUNKT , ERR );*)

   (*  MINDESTFEHLER DER INTERVALLMÄSSIGEN APPROXIMATION BERECHNEN  *)
   MINDIAM := SUP( ABS( CIFUNKT - MID( CIFUNKT ) ) );

   (*  JETZT TAYLOR-APPROXIMATION SUBTRAHIEREN                   *)
   (*  MITELWERTFORM, AUSGEWERTET NACH DEN REKURSIONSFORMEL AUS  *)
   (*  ACTON: "NUMERICAL METHODS THAT WORK"                      *)
   (*  VIC BERECHNET P( PHI_MID ), VID BERECHNET P'( PHI )       *)

   VIC[POLYORD+2] := 0;
   VIC[POLYORD+1] := 0;
   VID[POLYORD+2] := 0;
   VID[POLYORD+1] := 0;

   FOR I := POLYORD DOWNTO 0 DO BEGIN
       VIC[I] := ( 2 * COS( PHI_MID ) * VIC[I+1] - VIC[I+2] )
                   + ( RADIUS ** I ) * ##( TKOEFF[I].INT
                   + FOR L := 1 TO TKOEFF[I].PREC SUM( TKOEFF[I].STAG[L] ) );
       VID[I] := ( 2 * COS( PHI ) * VID[I+1] - VID[I+2] )
                   + ( RADIUS ** I ) * ##( I * TKOEFF[I].INT
                   + FOR L := 1 TO TKOEFF[I].PREC SUM( I*TKOEFF[I].STAG[L] ) );
   END;

   PINTMID := VIC[0] - VIC[1] * COS( PHI_MID )
                     + COMPL( 0 , 1 ) * VIC[1] * SIN( PHI_MID );

   PINT := PINTMID - ( PHI - PHI_MID ) * VID[1] * SIN( PHI )
           + COMPL(0,1) * ( PHI - PHI_MID ) * ( VID[0] - VID[1] * COS( PHI ) );

(* WRITELN( 'CIFUNKT, PINT' );
   WRITELN( CIFUNKT );
   WRITELN( PINT );
   WRITELN( 'CIMIDFUNKT, PINTMID' );
   WRITELN( CIMIDFUNKT );
   WRITELN( PINTMID ); *)

   CIFUNKT := CIFUNKT - PINT;
   CIMIDFUNKT := CIMIDFUNKT - PINTMID;

END;

(*----------------------------------------------------------------------*)

END.

(*==========================================================================*)
(*                                                                          *)
(*  PASCAL-XSC - MODUL  IVP_DATA                                            *)
(*                                                                          *)
(*  MODUL DER KOEFFIZIENTENFUNKTIONEN DES BEHANDELTEN ANFANGSWERTPROBLEMS   *)
(*                                                                          *)
(*==========================================================================*)
