(*==========================================================================*) (* *) (* PASCAL-XSC - MODUL PIVP_MPS 990602 *) (* *) (* MODUL ZUR L™SUNG DES ANFANGSWERTPROBLEMS ZUR DGL *) (* *) (* y^(n) = \sum_(i=0)^(n-1) p_i(x) y^(i) + p(x), p(x), p_i(x) Polynome *) (* *) (* POTENZREIHENANSATZ: y(x) = \sum_0^infty a_k x^k *) (* L™SUNG MIT VARIABLER TAYLORENTWICKLUNG UND RESTGLIEDABSCHŽTZUNG *) (* RECHNUNG MIT ADAPTIVER GENAUIGKEIT *) (* *) (* VERSION FšR RECORD-BASIERTE MPS-ARITHMETIK *) (* *) (*==========================================================================*) MODULE pivp_mps; USE iostd , i_ari , mv_ari , mvi_ari , service , mps_arip; GLOBAL TYPE INTEGERVECTOR = GLOBAL DYNAMIC ARRAY[*] OF INTEGER; INTEGERMATRIX = GLOBAL DYNAMIC ARRAY[*,*] OF INTEGER; (*--------------------------------------------------------------------------*) GLOBAL PROCEDURE IVP_POLY( N , M_AKT_MAX : INTEGER; VAR VRM : INTEGERVECTOR; VAR TIC : MPMATRIX; VAR VRY : RVECTOR; X_END , E_ABS : REAL; T_IND : INTEGER; VAR FAK : RVECTOR; VAR PRODINDEX : INTEGER; VAR PROD: RMATRIX; VAR IPROD: IMATRIX; VAR VIZ : IVECTOR; VAR OUTFILE , CONSOLE : TEXT ); (* N : ORNUNG DER DIFFERENTIALGLEICHUNG *) (* M_AKT_MAX : MAXIMALER GRAD DER POLYNOME *) CONST KMAX = 2000; (* MAXIMALE ORDNUNG DER TAYLORENTWICKLUNG VON Y *) (* VARIABLEN : *) (* AO_AKT : AKTUELLES MAXIMALES ZAHLENFORMAT *) (* PROD : HILSVARIABLE: PROD[I,J] := ( I+1 ) * ( I+2 ) *...* ( I+J ) *) (* TIA : TAYLORKOEFFIZIENTEN a_k * h^k DER LOESUNG *) (* TIC : TAYLORKOEFFIZIENTEN DER POLYNOME *) (* VIZ : FUNKTIONSWERTE DER N-1 ERSTEN ABLEITUNGEN AN DER STELLE X_END *) (* VRY : ANFANGSWERTE DER N-1 ERSTEN ABLEITUNGEN *) (* X_END : RECHTE GRENZE DES DEFINITIONSINTERVALLS DER DGL *) (* DO_KOEFSUM : ABLEITUNG VON KOEFSUM NACH OMEGA *) (* VORAUSSETZUNG AN PROD: DIE IN PROD BERECHNETEN PRODUKTE *) (* (MAXIMAL N + M_AKT_MAX FAKTOREN) MšSSEN EXAKT DARSTELLBAR SEIN *) (* FšR DIESE PRODUKTE K™NNTE MAN NOTFALLS EINE LANGZAHLDARSTELLUNG *) (* VERFšGBAR MACHEN *) VAR I , J , K , NR : INTEGER; AO_AKT : INTEGER; ANZAHL : INTEGER; AMAX , OMEGA : REAL; KOEFSUM , DO_KOEFSUM : REAL; RESTGLIED , DURCHMESSER : REAL; MAXSUMMAND , SKALFAKT : REAL; ABBRUCH : INTERVAL; UNTERLAUF , FINISHED : BOOLEAN; AUSDRUCKZEIGER : BOOLEAN; REIHE , IREST : IVECTOR[0..N]; MPX_END : MPVECTOR[0..N+M_AKT_MAX]; TIA : MPVECTOR[0..KMAX]; MINIM : INTEGERMATRIX[0..KMAX,0..N]; MULTIP : MPMATRIX[0..N,0..M_AKT_MAX]; IX_END : INTERVAL; (*--------------------------------------------------------------------------*) PROCEDURE AUSFILE1( K : INTEGER ); VAR I : INTEGER; BEGIN WRITELN( OUTFILE ); FOR I := 1 TO TIA[K].PREC DO BEGIN WRITE( OUTFILE , K : 3 , '. Summand = ' ); IF UNTERLAUF THEN BEGIN WRITE( OUTFILE , TIA[K].STAG[I] ); WRITELN( OUTFILE , ' (skal.)' ); END ELSE WRITELN( OUTFILE , TIA[K].STAG[I] ); END; IWRITELN( OUTFILE , TIA[K].INT ); WRITE( OUTFILE , K : 3 , '. Reihe = ' ); IF UNTERLAUF THEN BEGIN IWRITE( OUTFILE , REIHE[0] ); WRITELN( OUTFILE , ' (skal.)' ); END ELSE IWRITELN( OUTFILE , REIHE[0] ); WRITE( OUTFILE , K : 3 , '. Restglied : ' ); WRITELN( OUTFILE , RESTGLIED ); WRITE( OUTFILE , K : 3 , '. max. Wachstum : ' ); WRITELN( OUTFILE , KOEFSUM ); END; (*--------------------------------------------------------------------------*) PROCEDURE AUSFILE2( K : INTEGER ); VAR I : INTEGER; BEGIN WRITELN( OUTFILE ); FOR I := 1 TO TIA[K].PREC DO BEGIN WRITE( OUTFILE , K : 3 , '. Summand = ' ); IF UNTERLAUF THEN BEGIN WRITE( OUTFILE , TIA[K].STAG[I] ); WRITELN( OUTFILE , ' (skal.)' ); END ELSE WRITELN( OUTFILE , TIA[K].STAG[I] ); END; IWRITELN( OUTFILE , TIA[K].INT ); WRITELN( OUTFILE ); FOR I := 0 TO N - 1 DO BEGIN WRITE( OUTFILE , 'Restgl. fr y^(' , I : 1 , '): ' ); IF UNTERLAUF THEN BEGIN IWRITE( OUTFILE , IREST[I] ); WRITELN( OUTFILE , ' (skal.)' ); END ELSE IWRITELN( OUTFILE , IREST[I] ); END; FOR I := 0 TO N - 1 DO BEGIN WRITE( OUTFILE , 'Reihe fr y^(' , I : 1 , '): ' ); IF UNTERLAUF THEN BEGIN IWRITE( OUTFILE , REIHE[I] ); WRITELN( OUTFILE , ' (skal.)' ); END ELSE IWRITELN( OUTFILE , REIHE[I] ); END; END; (*--------------------------------------------------------------------------*) PROCEDURE AUSFILE3( K : INTEGER ); VAR I , J , NR : INTEGER; BEGIN WRITELN( OUTFILE ); WRITELN( OUTFILE, 'Zahlenformat nach ' , K:1 , ' Summanden erh”ht:' ); WRITELN( OUTFILE ); FOR NR := 0 TO K DO IF ( NR <= 23 ) OR ( ( NR <= 100 ) AND ( NR MOD 23 = 0 ) ) OR ( NR MOD 200 = 0 ) OR ( NR = K ) THEN BEGIN FOR I := 1 TO TIA[NR].PREC DO BEGIN WRITE( OUTFILE , NR : 3 , '. Summand = ' ); IF UNTERLAUF THEN BEGIN WRITE( OUTFILE , TIA[NR].STAG[I] ); WRITELN( OUTFILE , ' (skal.)' ); END ELSE WRITELN( OUTFILE , TIA[NR].STAG[I] / SKALFAKT ); END; IWRITELN( OUTFILE , TIA[NR].INT ); WRITELN( OUTFILE ); END; (* IF *) END; (*--------------------------------------------------------------------------*) PROCEDURE TAYKOEFF( K : INTEGER ); (* K+N-TEN TAYLORKOEFFIZIENTEN BERECHNEN *) VAR I , J , NR , PREC : INTEGER; DPL , DPU : DOTPRECISION; BEGIN PREC := TIA[K+N].PREC; IF K <= VRM[N] THEN BEGIN (* INHOMOGENITAET BERUECKSICHTIGEN *) DPL := #( FOR I := 1 TO TIC[N,K].PREC SUM( TIC[N,K].STAG[I] ) + TIC[N,K].INT.INF ); DPU := #( DPL + TIC[N,K].INT.SUP - TIC[N,K].INT.INF ); END ELSE BEGIN DPL := #( 0 ); DPU := DPL; END; FOR NR := 0 TO N - 1 DO FOR J := 0 TO MINIM[K,NR] DO ADD_MP_TIMES_MP( DPL , DPU , PROD[K-J,NR] *> MULTIP[NR,J] , TIA[K-J+NR] ); TIA[K+N] := CONVERT( DPL , DPU , PREC ); TIA[K+N].PREC := PREC; TIA[K+N] := TIA[K+N] / PROD[K,N]; WHILE ( ( PREC > 1 ) AND ( TIA[K+N].STAG[PREC] = 0 ) ) DO PREC := PREC - 1; TIA[K+N].PREC := PREC; END; (*--------------------------------------------------------------------------*) PROCEDURE FORMAT_CHECK; (* CHECK, OB DIE LŽNGE DES ZAHLENFORMATS ERH™HT WERDEN MUSS *) (* GLOBALE VARIABLEN : REIHE, N, K, TIA, VRY, PROD *) VAR I , J , J1 , J2 , NR , PREC : INTEGER; BEGIN IF ( ( REIHE[T_IND].SUP - REIHE[T_IND].INF ) / ( ABS(REIHE[T_IND].SUP ) + ABS( REIHE[T_IND].INF ) ) > 1.0E-14 ) AND ( REIHE[T_IND].SUP-REIHE[T_IND].INF > E_ABS ) THEN BEGIN (* ZAHLENFORMAT ERH™HEN *) (* MIT GENAUIGKEIT TIA[I].PREC+1 VON VORN BEGINNEN ! *) (* NULLTER SUMMAND IST MAXIMAL GENAU BERECHNET WORDEN *) (* SUMMANDEN 1 BIS (N-1) IM LANGZAHLFORMAT DARSTELLEN *) FOR I := 1 TO N-1 DO BEGIN PREC := TIA[I].PREC + 1; TIA[I] := ( MPX_END[I] *> VRY[I] ); IF TIA[I].PREC < PREC THEN TIA[I].PREC := PREC; TIA[I] := TIA[I] / FAK[I]; WHILE ( ( PREC > 1 ) AND ( TIA[I].STAG[PREC] = 0 ) ) DO PREC := PREC - 1; TIA[I].PREC := PREC; TIA[I].INT := TIA[I].INT * SKALFAKT; FOR J := 1 TO TIA[I].PREC DO TIA[I].STAG[J] := TIA[I].STAG[J] * SKALFAKT; END; (* WEITERE SUMMANDEN IM LANGZAHLFORMAT DARSTELLEN *) FOR I := 0 TO K DO BEGIN TIA[I+N].PREC := TIA[I+N].PREC + 1; TAYKOEFF( I ); (* FEHLERMELDUNG, FALLS UNTERLAUFBEREICH ERREICHT *) IF ( ( SUP( ABS( TIA[I+N].INT ) ) > 0 ) AND ( SUP( ABS( TIA[I+N].INT ) ) < 1.0E-308 ) ) THEN IF UNTERLAUF THEN BEGIN IF AUSDRUCKZEIGER THEN BEGIN WRITELN; WRITE( 'Unterlaufbereich bei dem ' , I + N : 1 ); WRITELN( '. Summanden erneut erreicht!'); WRITE( 'Eventuell keine maximal genaue Einschlieáung '); WRITELN( 'berechenbar!'); WRITELN; WRITELN( OUTFILE ); WRITE( OUTFILE , 'Unterlaufbereich bei dem ', I+N:1 ); WRITELN( OUTFILE , '. Summanden erneut erreicht!'); WRITE( OUTFILE , 'Eventuell keine maximal genaue '); WRITELN( OUTFILE , 'Einschlieáung berechenbar!'); WRITELN( OUTFILE ); AUSDRUCKZEIGER := FALSE; END ELSE BEGIN UNTERLAUF := TRUE; WRITELN; WRITE( 'Unterlaufbereich bei dem ' , I + N : 1 ); WRITELN( ' Summanden erreicht!'); WRITELN; WRITELN( OUTFILE ); WRITE( OUTFILE, 'Unterlaufbereich bei dem ', I + N : 1 ); WRITELN( OUTFILE , ' Summanden erreicht!'); WRITELN( OUTFILE ); FOR J := 1 TO TIA[I+N].PREC DO BEGIN WRITE( OUTFILE , I + N : 3 , '. Summand = ' ); IF UNTERLAUF THEN BEGIN WRITE( OUTFILE , TIA[I+N].STAG[J] ); WRITELN( OUTFILE , ' (skal.)' ); END ELSE BEGIN WRITELN( OUTFILE , TIA[I+N].STAG[J] ); END; END; (* FOR J *) (* SKALIEREN, UM DEN GANZEN DARSTELLBAREN *) (* ZAHLENBEREICH DES RECHNERS ZU NUTZEN *) IF KOEFSUM < 1 THEN BEGIN (* SKALIERUNG NUR DURCHFUEHREN, *) (* FALLS WACHSTUMSPHASE BEENDET *) ANZAHL := TRUNC( LOG10( 1.0E+307/MAXSUMMAND ) / 9 ) - 1; SKALFAKT := 1073741824; FOR J1 := 1 TO ANZAHL DO SKALFAKT := SKALFAKT * 1073741824; WRITELN; WRITELN( 'Skalierung mit Faktor ' , SKALFAKT ); WRITELN; WRITELN( OUTFILE ); WRITE( OUTFILE , 'Skalierung mit Faktor ', SKALFAKT ); WRITELN( OUTFILE ); FOR J1 := 0 TO K + N DO BEGIN TIA[J1].INT := TIA[J1].INT * SKALFAKT; FOR J2 := 1 TO TIA[J1].PREC DO TIA[J1].STAG[J2] := TIA[J1].STAG[J2] * SKALFAKT; END; (* INHOMOGENITAET SKALIEREN *) FOR J1 := VRM[N] DOWNTO 0 DO BEGIN TIC[N,J1].INT := TIC[N,J1].INT * SKALFAKT; FOR J2 := TIC[N,J1].PREC DOWNTO 0 DO TIC[N,J1].STAG[J2]:=TIC[N,J1].STAG[J2]*SKALFAKT; END; END; (* IF KOEFSUM < 1 *) END; (* IF AUSDRUCKZEIGER *) END; (* IF UNTERLAUF *) END; (* FOR I := 0 TO K *) (* WERTE AUSGEBEN *) {REIHE[T_IND] := ##( FOR I := T_IND TO K + N SUM( FOR J := 1 TO TIA[I].PREC SUM( PROD[I-T_IND,T_IND] * TIA[I].STAG[J] ) ) ); REIHE[T_IND] := REIHE[T_IND] / ( ##( FOR NR := 1 TO MPX_END[T_IND].PREC SUM( MPX_END[T_IND].STAG[NR] ) + MPX_END[T_IND].INT ) ); WRITELN( 'Zahlenformat nach ' , K : 1 , ' Summanden erh”ht.' ); AUSFILE3( K + N );} (* AKTUELL BENOETIGTE GENAUIGKEIT ERGIBT SICH *) (* AUS DER GENAUIGKEIT DER LETZTEN N SUMMANDEN *) AO_AKT := TIA[K+1].PREC; FOR I := 2 TO N DO IF TIA[K+I].PREC > AO_AKT THEN AO_AKT := TIA[K+I].PREC; END; (* IF DIAM( REIHE ) ZU GROSS *) END; (* FORMAT_CHECK *) (*--------------------------------------------------------------------------*) PROCEDURE PRODUCTS( K : INTEGER; VAR PRODINDEX : INTEGER ); (* (K+N)-TE PRODUKTE BERECHNEN *) VAR I , J : INTEGER; BEGIN FOR I := 0 TO N DO BEGIN IPROD[K+N-I,I] := IPROD[K+N-I-1,I] * ( K + N ) / ( K + N - I ); IF IPROD[K+N-I,I].INF < IPROD[K+N-I,I].SUP THEN BEGIN WRITELN; WRITELN( 'Fehler in Modul PIVP_MPS' ); WRITELN( 'IPROD in Prozedur PRODUCTS nicht exakt darstellbar.' ); WRITELN; END ELSE PROD[K+N-I,I] := IPROD[K+N-I,I].INF; END; PRODINDEX := K; END; (* *) (*--------------------------------------------------------------------------*) (* *) (* HAUPTPROGRAM *) (* *) (*--------------------------------------------------------------------------*) (* *) BEGIN (* ZUERST POTENZEN VON X_END BEREITSTELLEN *) MPX_END[0].PREC := 1; MPX_END[0].STAG := 1; MPX_END[0].INT := 0; MPX_END[1].PREC := 1; MPX_END[1].STAG := X_END; MPX_END[1].INT := 0; FOR I := 2 TO N + M_AKT_MAX DO MPX_END[I] := MPX_END[I-1] *> X_END; (* STARTWERTE VORBESETZEN *) UNTERLAUF := FALSE; AUSDRUCKZEIGER := TRUE; SKALFAKT := 1; TIA[0].PREC := 1; TIA[0].STAG := VRY[0]; TIA[0].INT := 0; AO_AKT := 1; FOR I := 1 TO N - 1 DO BEGIN TIA[I] := ( MPX_END[I] *> VRY[I] ) / FAK[I]; TIA[I] := ROUND( TIA[I] , 1 ); END; (* BASIS-MULTIPLIKATOREN BERECHNEN *) FOR I := 0 TO N - 1 DO FOR J := 0 TO VRM[I] DO MULTIP[I,J] := TIC[I,J] *> MPX_END[N-I+J]; K := - 1; KOEFSUM := 1; WHILE ( K <= N + M_AKT_MAX ) OR ( ( KOEFSUM >= 0.97 ) AND ( K < KMAX - N ) ) DO BEGIN (* "VORITERATION" OHNE FEHLERABSCHŽTZUNG *) (* KOEFSUM GIBT AN, OB DIE a_k NOCH WACHSEN K™NNEN *) (* GETESTET WIRD FšR OMEGA = 1, UM RECHENZEIT ZU SPAREN *) (* ABBRUCH DAFšR ERST, WENN KOEFSUM < 0.97 STATT 1.00 *) K := K + 1; TIA[K+N].PREC := AO_AKT; (* K+N-TEN TAYLORKOEFFIZIENTEN BERECHNEN *) FOR I := 0 TO N-1 DO MINIM[K,I] := MIN( K , VRM[I] ); TAYKOEFF( K ); (* (K+N+1)-TE PRODUKTE BERECHNEN *) IF K + 1 > PRODINDEX THEN PRODUCTS( K + 1 , PRODINDEX ); (* KONTRAKTIONSBEDINGUNG PRšFEN *) (* šBERPRšFT WIRD DIE KONTRAKTIONSBEDINGUNG *) (* FšR DIE N-1-te ABLEITUNG *) (* DURCH STAGGERED CORRECTION IST DIE SUMME DER *) (* MULTIP[K,I,J].STAG[L] FšR L > 0 BETRAGSMŽSSIG *) (* KLEINER ALS MULTIP[K,I,J].STAG[1] * 10^-12 *) (* UM RECHENZEIT ZU SPAREN, WIRD KEINE GENAUERE *) (* PRšFUNG DURCHGEFšHRT, SONDERN KOEFSUM STETS *) (* MIT ( 1 + 10^-12 ) MULTIPLIZIERT *) (* (GERINGE šBERSCHŽTZUNG OHNE PRAKTISCH *) (* RELEVANTEN EINFLUSS AUF DIE ABBRUCHBEDINGUNG) *) IF K > N + M_AKT_MAX THEN BEGIN (* KOEFSUM WIRD FšR DEN INDEX K+N+1 BERECHNET *) KOEFSUM := 0; FOR I := 0 TO N - 1 DO FOR J := 0 TO VRM[I] DO KOEFSUM := KOEFSUM + PROD[K+1-J,I] * ( ABS( MULTIP[I,J].STAG[1] ) + SUP( ABS( MULTIP[I,J].INT ) ) ); (* KONTRAKTIONSBEDINGUNG FšR DIE N-1-te ABLEITUNG *) (* ERFORDERT DIVISION DURCH PROD[K+1-M_AKT_MAX-(N-1),N] *) (* ANSTELLE VON PROD[K+1,N] *) KOEFSUM := 1.000000000001 * KOEFSUM / PROD[K-M_AKT_MAX-N+2,N]; END; (* CHECK, OB DIE LŽNGE DES ZAHLENFORMATS ERH™HT WERDEN MUSS *) (* ZU HŽUFIGEN TEST VERHINDERN *) IF ( K > N ) AND ( K MOD 3 = 0 ) THEN BEGIN REIHE[T_IND] := ##( FOR I := T_IND TO K + N SUM( PROD[I-T_IND,T_IND] * TIA[I].INT + FOR J := 1 TO TIA[I].PREC SUM( PROD[I-T_IND,T_IND] * TIA[I].STAG[J] ))); REIHE[T_IND] := REIHE[T_IND] / ( ##( FOR NR := 1 TO MPX_END[T_IND].PREC SUM( MPX_END[T_IND].STAG[NR] )+MPX_END[T_IND].INT ) ); FORMAT_CHECK; END; {IF ( K + N <= 23 ) OR ( ( K + N <=100 ) AND ( ( K+N ) MOD 23 = 0 ) ) OR ( ( K + N ) MOD 200 = 0 ) THEN BEGIN REIHE[0] := ##( FOR I := 0 TO K + N SUM( FOR J := 1 TO TIA[I].PREC SUM( TIA[I].STAG[J] ) ) ); AUSFILE1( K + N ); END;} (* IF ( K + N ) MOD 89 = 0 THEN IF UNTERLAUF THEN WRITELN( CONSOLE , K+N : 3 , ' ' , TIA[K+N].STAG[1] , ' (skal.)' ) ELSE WRITELN( CONSOLE , K + N : 3 , ' ' , TIA[K+N].STAG[1] ); FLUSH( CONSOLE );*) END; (* WHILE *) WRITELN; WRITE( 'Wachstum der Summanden nach ' , K + N : 1 ); WRITELN( ' Summanden beendet.' ); {WRITELN( OUTFILE , ); WRITE( OUTFILE , 'Wachstum der Summanden nach ' , K + N : 1 ); WRITELN( OUTFILE , ' Summanden beendet.' );} (* MAXSUMMAND: NUR ZUR SKALIERUNG BEI NUMERISCHEM šBERLAUF *) MAXSUMMAND := ABS( TIA[K+N].STAG[1] ); (* ITERATION MIT FEHLERABSCHŽTZUNG *) (* HIER IST K SO GROSS, DASS DIE *) (* INHOMOGENITŽT NICHT MEHR AUFTRITT *) (* HILFPARAMETER OMEGA "HEURISTISCH OPTIMAL" WŽHLEN *) (* FšR OMEGA = 1 GILT MOMENTAN KOEFSUM = 0.97 *) (* NEWTON-VERBESSERUNG VON OMEGA, SO DASS KOEFSUM = 0.98 *) DO_KOEFSUM := 0; FOR I := 0 TO N - 1 DO FOR J := 0 TO VRM[I] DO DO_KOEFSUM := DO_KOEFSUM - PROD[K+1-J,I] * ( N - I + J ) * ( ABS( MULTIP[I,J].STAG[1] ) + SUP( ABS( MULTIP[I,J].INT ) ) ); (* KONTRAKTIONSBEDINGUNG FšR DIE N-1-te ABLEITUNG *) (* ERFORDERT DIVISION DURCH PROD[K+1-M_AKT_MAX-(N-1),N] *) (* ANSTELLE VON PROD[K+1,N] *) DO_KOEFSUM := 1.000000000001 * DO_KOEFSUM / PROD[K-M_AKT_MAX-N+2,N]; OMEGA := 1 - ( KOEFSUM - 0.98 ) / DO_KOEFSUM; {WRITELN( OUTFILE , 'OMAGA = ' , OMEGA ); WRITELN( OUTFILE , 'DO_KOEFSUM = ' , DO_KOEFSUM );} (* TESTEN: OMEGA \IN (0,1)? *) IF ( OMEGA < 0 ) OR ( OMEGA >= 1 ) THEN BEGIN WRITELN; WRITELN( 'Fehler: OMAGA ist nicht geignet!' ); WRITELN( 'OMAGA = ' , OMEGA ); WRITELN( 'DO_KOEFSUM = ' , DO_KOEFSUM ); (* ABBRUCH ERZEUGEN *) ABBRUCH := INTVAL( 1.23 , -1.23 ); END; REPEAT K := K + 1; TIA[K+N].PREC := AO_AKT; (* (K+N)-TEN TAYLORKOEFFIZIENTEN BERECHNEN *) FOR I := 0 TO N-1 DO MINIM[K,I] := VRM[I]; (* = MIN( K , VRM[I] ) *) TAYKOEFF( K ); (* (K+N+1)-TE PRODUKTE BERECHNEN *) IF K + 1 > PRODINDEX THEN PRODUCTS( K + 1 , PRODINDEX ); (* IF ( K + N ) MOD 89 = 0 THEN IF UNTERLAUF THEN WRITELN( CONSOLE, K+N:3, ' ', TIA[K+N].STAG[1], ' (skal.)' ) ELSE WRITELN( CONSOLE, K+N:3, ' ', TIA[K+N].STAG[1] ); FLUSH( CONSOLE ); *) (* ABBRUCHKRITERIUM: *) (* T_IND-te ABLEITUNG MAXIMAL GENAU EINGESCHLOSSEN *) FINISHED := TRUE; FOR I := T_IND TO T_IND DO BEGIN (* I-TE ABLEITUNG EINSCHLIESSEN *) REIHE[I] := ##( FOR J := I TO K + N SUM( PROD[J-I,I] * TIA[J].INT + FOR NR := 1 TO TIA[J].PREC SUM( PROD[J-I,I] * TIA[J].STAG[NR] ) ) ); (* RESTGLIED ABSCHŽTZEN *) (* ZUERST MAXIMUM DER LETZTEN N KOEFFIZIENTEN BESTIMMEN *) (* DURCH STAGGERED CORRECTION IST DIE SUMME DER *) (* TIA[K].STAG[I] FšR I > 0 BETRAGSMŽSSIG KLEINER ALS *) (* TIA[K].STAG[1] * 10^-12 *) (* UM RECHENZEIT ZU SPAREN, WIRD KEINE GENAUERE PRšFUNG *) (* DURCHGEFšHRT, SONDERN AMAX STETS MIT ( 1 + 10^-12 ) *) (* MULTIPLIZIERT (GERINGE šBERSCHŽTZUNG DES RESTGLIEDS OHNE *) (* PRAKTISCH RELEVANTEN EINFLUSS AUF DIE ABBRUCHBEDINGUNG) *) AMAX := ABS( PROD[K+N-I,I] * ( TIA[K+N].STAG[1] + SUP( ABS( TIA[K+N].INT ) ) ) ); FOR J := 1 TO MIN( K , N + M_AKT_MAX - 1 ) DO IF ABS( PROD[K+N-J-I,I] * ( TIA[K+N-J].STAG[1] + SUP( ABS( TIA[K+N-J].INT ))) * ( OMEGA ** J ) ) > AMAX THEN AMAX := ABS( PROD[K+N-J-I,I] * ( TIA[K+N-J].STAG[1] + SUP( ABS( TIA[K+N-J].INT ))) * ( OMEGA ** J ) ); AMAX := AMAX * 1.000000000001; (* KOEFSUM FšR DEN INDEX K+N+1 BERECHNEN *) KOEFSUM := 0; FOR NR := 0 TO N - 1 DO FOR J := 0 TO VRM[NR] DO KOEFSUM := KOEFSUM + PROD[K+1-J,NR] * ( ABS( MULTIP[NR,J].STAG[1] ) + SUP( ABS( MULTIP[NR,J].INT ) ) ) / ( OMEGA ** (N-NR+J) ); (* KONTRAKTIONSBEDINGUNG FšR DIE N-1-te ABLEITUNG *) (* ERFORDERT DIVISION DURCH PROD[K+1-M_AKT_MAX-(N-1),N] *) (* ANSTELLE VON PROD[K+1,N] *) KOEFSUM := 1.000000000001 * KOEFSUM / PROD[K-M_AKT_MAX-N+2,N]; {WRITELN( OUTFILE , 'KOEFSUM = ' , KOEFSUM );} (* PRšFEN: KOEFSUM \IN (0,1) ? *) IF ( KOEFSUM < 0 ) OR ( KOEFSUM >= 1 ) THEN BEGIN WRITELN( 'Fehler: KOEFSUM ist nicht geignet!' ); WRITELN( 'KOEFSUM = ' , KOEFSUM ); (* ABBRUCH ERZEUGEN *) ABBRUCH := INTVAL( 1.23 , -1.23 ); END; RESTGLIED := AMAX * OMEGA / ( 1 - OMEGA ); IREST[I] := INTVAL( - RESTGLIED , RESTGLIED ); IX_END := ##( FOR NR := 1 TO MPX_END[I].PREC SUM( MPX_END[I].STAG[NR] ) + MPX_END[I].INT ); VIZ[I] := ( REIHE[I] + IREST[I] ) / IX_END; REIHE[I] := REIHE[I] / IX_END; IREST[I] := IREST[I] / IX_END; DURCHMESSER := VIZ[I].SUP -> VIZ[I].INF; IF ( RESTGLIED / ABS( REIHE[I].SUP ) > 1.0E-17 ) AND ( DURCHMESSER > E_ABS ) THEN FINISHED := FALSE; {AUSFILE1( K + N );} END; IF NOT FINISHED THEN BEGIN FORMAT_CHECK; (* NEWTON-VERBESSERUNG VON OMEGA, SO DASS KOEFSUM = 0.98 *) DO_KOEFSUM := 0; FOR I := 0 TO N - 1 DO FOR J := 0 TO VRM[I] DO DO_KOEFSUM := DO_KOEFSUM - PROD[K+1-J,I] * (N-I+J) * ( ABS( MULTIP[I,J].STAG[1] ) + SUP( ABS( MULTIP[I,J].INT ) ) ) / ( OMEGA ** (N-I+J+1) ); (* KONTRAKTIONSBEDINGUNG FšR DIE N-1-te ABLEITUNG *) (* ERFORDERT DIVISION DURCH PROD[K+1-M_AKT_MAX-(N-1),N] *) (* ANSTELLE VON PROD[K+1,N] *) DO_KOEFSUM := 1.000000000001*DO_KOEFSUM / PROD[K-M_AKT_MAX-N+2,N]; OMEGA := OMEGA - ( KOEFSUM - 0.98 ) / DO_KOEFSUM; {WRITELN( OUTFILE , 'DO_KOEFSUM = ' , DO_KOEFSUM ); WRITELN( OUTFILE , 'OMAGA = ' , OMEGA );} (* TESTEN: OMEGA \IN (0,1)? *) IF ( OMEGA < 0 ) OR ( OMEGA >= 1 ) THEN BEGIN WRITELN( 'Fehler: OMAGA ist nicht geignet!' ); WRITELN( 'DO_KOEFSUM = ' , DO_KOEFSUM ); WRITELN( 'OMAGA = ' , OMEGA ); (* ABBRUCH ERZEUGEN *) ABBRUCH := INTVAL( 1.23 , -1.23 ); END; END; UNTIL FINISHED; FOR I := N - 1 DOWNTO 0 DO BEGIN (* I-TE ABLEITUNG EINSCHLIESSEN *) REIHE[I] := ##( FOR J := I TO K + N SUM( PROD[J-I,I] * TIA[J].INT + FOR NR := 1 TO TIA[J].PREC SUM( PROD[J-I,I] * TIA[J].STAG[NR] ) ) ); (* RESTGLIED ABSCHŽTZEN *) (* ZUERST MAXIMUM DER LETZTEN N KOEFFIZIENTEN BESTIMMEN *) (* DURCH STAGGERED CORRECTION IST DIE SUMME DER TIA[K] *) (* FšR I > 0 BETRAGSMŽSSIG KLEINER ALS *) (* TIA[K].STAG[1] * 10^-12 *) (* UM RECHENZEIT ZU SPAREN, WIRD KEINE GENAUERE PRšFUNG *) (* DURCHGEFšHRT, SONDERN AMAX STETS MIT ( 1 + 10^-12 ) *) (* MULTIPLIZIERT (GERINGE šBERSCHŽTZUNG DES RESTGLIEDS OHNE *) (* PRAKTISCH RELEVANTEN EINFLUSS AUF DIE ABBRUCHBEDINGUNG) *) AMAX := ABS( PROD[K+N-I,I] * ( TIA[K+N].STAG[1] + SUP( ABS( TIA[K+N].INT ) ) ) ); FOR J := 1 TO MIN( K , N + M_AKT_MAX - 1 ) DO IF ABS( PROD[K+N-J-I,I] * ( TIA[K+N-J].STAG[1] + SUP( ABS( TIA[K+N-J].INT ))) * ( OMEGA ** J ) ) > AMAX THEN AMAX := ABS( PROD[K+N-J-I,I] * ( TIA[K+N-J].STAG[1] + SUP( ABS( TIA[K+N-J].INT ))) * ( OMEGA ** J ) ); AMAX := AMAX * 1.000000000001; RESTGLIED := AMAX * OMEGA / ( 1 - OMEGA ); IREST[I] := INTVAL( - RESTGLIED , RESTGLIED ); IX_END := ##( FOR NR := 1 TO MPX_END[I].PREC SUM( MPX_END[I].STAG[NR] ) + MPX_END[I].INT ); VIZ[I] := ( REIHE[I] + IREST[I] ) / IX_END; REIHE[I] := REIHE[I] / IX_END; IREST[I] := IREST[I] / IX_END; END; IF SKALFAKT <> 1 THEN BEGIN (* RUECKSKALIERUNG *) FOR I := N - 1 DOWNTO 0 DO VIZ[I] := VIZ[I] / SKALFAKT; (* INHOMOGENITAET RUEKSKALIEREN *) FOR J := VRM[N] DOWNTO 0 DO BEGIN TIC[N,J].INT := TIC[N,J].INT / SKALFAKT; FOR NR := TIC[N,J].PREC DOWNTO 0 DO TIC[N,J].STAG[NR] := TIC[N,J].STAG[NR] / SKALFAKT; END; END; WRITELN( 'Einschlieáung mit ' , K + N : 1 , ' Reihengliedern erhalten.' ); {AUSFILE2( K + N ); WRITELN( OUTFILE ); WRITELN( OUTFILE ); WRITELN( OUTFILE , 'Einschlieáung mit ' , K + N : 1 , ' Reihengliedern:' ); WRITELN( OUTFILE ); FOR I := 0 TO N - 1 DO BEGIN WRITE( OUTFILE , 'y^(' , I : 1 , ') (' , X_END : 10 : 9 , ') î ' ); IWRITELN( OUTFILE , VIZ[I] ); END;} END; END. (*==========================================================================*) (* *) (* PASCAL-XSC - MODUL PIVP_MPS *) (* *) (* MODUL ZUR L™SUNG DES ANFANGSWERTPROBLEMS ZUR DGL *) (* *) (* y^(n) = \sum_(i=0)^(n-1) p_i(x) y^(i) + p(x), p(x), p_i(x) Polynome *) (* *) (* POTENZREIHENANSATZ: y(x) = \sum_0^infty a_k x^k *) (* L™SUNG MIT VARIABLER TAYLORENTWICKLUNG UND RESTGLIEDABSCHŽTZUNG *) (* RECHNUNG MIT ADAPTIVER GENAUIGKEIT *) (* *) (* VERSION FšR RECORD-BASIERTE MPS-ARITHMETIK *) (* *) (*==========================================================================*)