(*==========================================================================*) (* *) (* PASCAL-XSC - PROGRAMM LIVPTAYP 000222 *) (* *) (* PROGRAMM 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) analytisch *) (* *) (* POTENZREIHENANSATZ: y(x) = \sum_0^infty a_k x^k *) (* LÖSUNG MIT VARIABLER TAYLORENTWICKLUNG UND RESTGLIEDABSCHÄTZUNG *) (* RECHNUNG MIT ADAPTIVER GENAUIGKEIT *) (* *) (* HIER : HAUPTPROGRAMM *) (* EINLESEN DES ANFANGSWERTPROBLEMS *) (* ZERLEGEN DES DEFINITIONSINTERVALLS IN GEEIGNETE ABSCHNITTE *) (* AUFRUF DER PROZEDUR LIVP IM MODUL LIVP_MPS *) (* *) (*==========================================================================*) PROGRAM LIVPTAYP( INPUT , OUTPUT ); USE iostd , x_real, i_ari , mv_ari , mvi_ari , service , timer , lss , ilss , mps_aril , itaylor , mps_tayl , ivp_data , cauchyco , livp_mps; CONST KMAX = 300 ; (* MAXIMALER GRAD DER TAYLORENTWICKLUNG *) (* T_IND : TEST-INDEX FÜR ABBRUCHKRITERIUM IM UNTERPROGRAMM LIVP *) (* 0 FÜR FUNKTION SELBST, J FÜR J-TE ABLEITUNG *) (* ORDNUNG : ORDNUNG DER POLYNOM-ANTEILE DER TAYLORREIHEN *) VAR N , BSPNR : INTEGER; CPU_TIME : REAL; AKT_DATUM , AWP_STRING , LOESUNG : STRING; DATEIN , OUTFILE , CONSOLE : TEXT; (*--------------------------------------------------------------------------*) PROCEDURE EINLESEN_MAIN; BEGIN RESET( INPUT , 'awp.ein' ); REWRITE( OUTPUT , 'awp.aus' ); REWRITE( OUTFILE , 'outfile.awp' ); REWRITE( CONSOLE , 'con' ); (* TAGESDATUM EINLESEN AUS FILE *) RESET( DATEIN , 'aktdatum.txt' ); READLN( DATEIN , AKT_DATUM ); (* FUNKTIONSNAME EINLESEN *) READLN; READLN; READLN( AWP_STRING ); READLN; READLN( LOESUNG ); READLN; (* NUMMER DES BEISPIELS EINLESEN *) READLN( BSPNR ); (* ORDNUNG DER DGL EINLESEN *) READLN( N ); END; (*--------------------------------------------------------------------------*) PROCEDURE IVP_HAUPT( N : INTEGER ); (* VARIABLEN : *) (* TIB : TAYLORKOEFFIZIENTEN DER KOEFFIZIENTENFUNKTIONEN *) (* VRM : ORDNUNG DER APPROXIMATIONSPOLYNOME DER KOEFFIZIENTENFUNKTIONEN *) (* VRB : CAUCHY-SCHRANKE DER TAYLORKOEFFIZIENTEN MIT ORDNUNG > VRM *) (* X_END : RECHTE GRENZE DES DEFINITIONSINTERVALLS DER DGL *) VAR I , J , K : INTEGER; KAPPA , SCHRITTE , T_IND , VARIANTEN : INTEGER; ERR , SUBINTERVALS , KOEFPREC : INTEGER; ORDNUNG , INHDGL : INTEGER; SEKTOREN , SEKTOR_MAX : INTEGER; INHOM : BOOLEAN; M_AKT_MAX , AW_REAL , INH_GRAD : INTEGER; ALPHA , BETA , GAMMA , OMEGA : REAL; E_ABS , RMAX_MID_LB , RMAX_LB : REAL; RMAX , X_AKT , X_END : REAL; VRM : INTEGERVECTOR[-1..N-1]; VRY , VRNULL : RVECTOR[0..N-1]; VR_PART1 , VR_PART4 , VR_PART5 : RVECTOR[0..N-1]; VIY , VIZ , VI_INH , VINULL : IVECTOR[0..N-1]; VI_PART1 , VI_PART4 , VI_PART5 : IVECTOR[0..N-1]; VIY0 , VIY1 : IVECTOR[0..N-1]; VIY2 , VIY3 , VIY4 , VIY5 : IVECTOR[0..N-1]; VIU2 , VIU3 , VIU4 , VIU5 : IVECTOR[0..N-1]; MRID : RMATRIX[0..N-1,0..N-1]; MRB , MRB2 , MRB3 , MRB4 , MRB5 : RMATRIX[0..N-1,0..N-1]; TIB : MPMATRIX[-1..N-1,0..KMAX]; VRB : RVECTOR[-1..N-1]; MIA , MIBINV2 , MIBINV3 : IMATRIX[0..N-1,0..N-1]; MIBINV4 , MIBINV5 : IMATRIX[0..N-1,0..N-1]; PRODINDEX : INTEGER; FAK : RVECTOR[0..30]; PROD : RMATRIX[0..KMAX,0..N]; IPROD : IMATRIX[0..KMAX,0..N]; (*--------------------------------------------------------------------------*) PROCEDURE EINLESEN; VAR I , J : INTEGER; BEGIN (* HOMOGENE ODER INHOMOGENE DGL? *) READLN( INHDGL ); IF INHDGL = 0 THEN INHOM := FALSE ELSE INHOM := TRUE; (* PARAMETER DER DGL *) READLN( ALPHA ); READLN( BETA ); READLN( GAMMA ); (* ANFANGSWERTE EINLESEN *) READLN( AW_REAL ); IF AW_REAL = 1 THEN (* REELLE ANFANGSWERTE *) FOR I := 0 TO N - 1 DO BEGIN READLN( VRY[I] ); VIY[I] := VRY[I]; END ELSE IF AW_REAL = 2 THEN (* REELLE HEXADEZIMALE ANFANGSWERTE *) FOR I := 0 TO N - 1 DO BEGIN READLN( VRY[I] : 'X' ); VIY[I] := VRY[I]; END ELSE IF AW_REAL = 3 THEN (* INTERVALL-ANFANGSWERTE *) FOR I := 0 TO N - 1 DO READLN( VIY[I] ) ELSE (* HEXADEZIMALE INTERVALL-ANFANGSWERTE *) FOR I := 0 TO N - 1 DO BEGIN READ( VRY[I] : 'X' ); VIY[I].INF := VRY[I]; READLN( VRY[I] : 'X' ); VIY[I].SUP := VRY[I]; END; (* RECHTE GRENZE DES DEFINITIONSINTERVALLS DER DGL EINLESEN *) READLN( X_END ); (* GEWÜNSCHTE GENAUIGKEIT DER TAYLORKOEFFIZIENTEN *) (* DER KOEFFIZIENTENFUNKTIONEN EINLESEN UND *) (* GENAUIGKEITSVARIABLE FÜR TAYLOR-MODUL SETZEN *) READLN( KOEFPREC ); SET_TAYPREC( KOEFPREC ); SETPREC( 3 * KOEFPREC + 3 ); (* MAXIMALE ANZAHL DER SEKTOREN FÜR DIE *) (* CAUCHYSCHE KOEFFIZIENTENABSCHÄTZUNG EINLESEN *) READLN( SEKTOR_MAX ); (* OMEGA FÜR DIE GEOMETRISCHE REIHE EINLESEN *) READLN( OMEGA ); (* ABSOLUTE FEHLERSCHRANKE EINLESEN *) (* NUR FÜR UNTERLAUFBEREICH INTERESSANT *) READLN( E_ABS ); (* TEST-INDEX EINLESEN *) (* T_IND-te ABLEITUNG WIRD AUF GENAUIGKEIT ÜBERPRÜFT *) READLN( T_IND ); (* SCHRITTE EINLESEN (1 FÜR EINEN SCHRITT, 0 FÜR GESTEUERTE SCHRITTWEITE *) READLN( SCHRITTE ); (* KAPPA FÜR REZESSIONSBEDINGUNG EINLESEN *) (* FALLS DIE REZESSIONSBEDINGUNG NACH KAPPA SUMMANDEN *) (* NICHT ERFÜLLT IST, WIRD DAS INTERVALL HALBIERT *) (* NUR BEI GESTEUERTER SCHRITTWEITE BENUTZT *) READLN( KAPPA ); (* EINSCHLIESSUNGSVARIANTEN EINLESEN *) READLN( VARIANTEN ); END; (*--------------------------------------------------------------------------*) PROCEDURE WRITE_HEAD; VAR I : INTEGER; BEGIN WRITE( 'Verifizierte Lösung des Anfangswertproblems ' ); WRITELN( ' ' , AKT_DATUM ); WRITELN; WRITELN( AWP_STRING ); WRITELN; FOR I := 0 TO N - 1 DO IF VIY[I].INF = VIY[I].SUP THEN BEGIN WRITE( 'y^(' , I : 1 , ') (0) = ' , VIY[I].SUP , ' = ' ); WRITELN( VIY[I].SUP : 'X' ); END ELSE BEGIN WRITE( 'y^(' , I : 1 , ') (0) î ' ); IWRITELN( VIY[I] ); WRITE( ' = ' ); WRITELN( '[ ' , VIY[I].INF : 'X', ' , ', VIY[I].SUP : 'X' , ' ]' ); END; WRITELN; WRITELN( 'im Intervall [ 0 , ' , X_END , ' ]' ); WRITELN; WRITELN( LOESUNG ); IF BSPNR = 1 THEN BEGIN WRITELN; WRITELN( 'ALPHA = ' , ALPHA ); END; IF ( BSPNR = 2 ) OR ( BSPNR = 4 ) THEN BEGIN WRITELN; WRITELN( 'ALPHA = ', ALPHA , ' BETA = ' , BETA , ' GAMMA = ' , GAMMA ); END; WRITELN; WRITELN( 'Programm LIVPTAYP vom 22.02.00' ); WRITELN; WRITELN( 'Potenzreihenansatz mit adaptiver Genauigkeit.' ); IF SCHRITTE = 1 THEN WRITELN( 'Rechnung mit einem Integrationsschritt.' ) ELSE BEGIN WRITE( 'Schrittweitensteuerung mit Kontraktionsbedingung nach ' ); WRITELN( KAPPA : 1 , ' Koeffizienten.' ); END; WRITELN( 'Rechengenauigkeit: ' , KOEFPREC + 1 : 1 , '.' ); WRITELN( 'Anzahl Sektoren : ' , SEKTOR_MAX , '.' ); WRITELN; WRITELN( 'Abbruchkriterium: omega = ' , OMEGA , '.' ); WRITELN( ' E_abs = ' , E_ABS , '.' ); WRITE( ' Überprüft wurde die Genauigkeit der ' ); WRITELN( T_IND : 1 , '. Ableitung.' ); WRITELN; END; (*--------------------------------------------------------------------------*) PROCEDURE QR( VAR MRA : RMATRIX ); (* QR-ZERLEGUNG DER MATRIX MRA *) (* AUSGABE IM PARAMETER MRA : ORTHOGONALE MATRIX MRQ *) VAR I , J , K , NL , NU : INTEGER; MRQ : RMATRIX[LBOUND(MRA)..UBOUND(MRA),LBOUND(MRA)..UBOUND(MRA)]; B,C : RVECTOR[LBOUND(MRA)..UBOUND(MRA)]; S : REAL; BEGIN NL := LBOUND(MRA); NU := UBOUND(MRA); MRQ := ID( MRQ ); FOR K := NL TO NU-1 DO BEGIN S := 0; FOR I := K TO NU DO BEGIN B[I] := MRA[I,K]; IF S < ABS( B[I] ) THEN S := ABS( B[I] ); END; IF S > 0 THEN BEGIN S := SQRT( B * B ); IF MRA[K,K] >= 0 THEN S := -S; B[K] := MRA[K,K] - S; S := 1 / ( S * B[K] ); C := S * B; FOR J := K + 1 TO NU DO BEGIN S := C * RVECTOR( MRA[*,J] ); FOR I := K TO NU DO MRA[I,J] := MRA[I,J] + B[I]*S; END; FOR I := NL TO NU DO BEGIN S := RVECTOR( MRQ[I] ) * C; FOR J := K TO NU DO MRQ[I,J] := MRQ[I,J] + S * B[J]; END; END; B[K] := 0; C[K] := 0; END; MRA := MRQ; END; (*--------------------------------------------------------------------------*) PROCEDURE QR_SORT( VAR MRA : RMATRIX; REST : IVECTOR ); (* QR-ZERLEGUNG VON MRA MIT SORTIEREN DER SPALTEN VON MRA *) (* GLOBALE VARIABLE : N *) VAR I , J , K : INTEGER; HV , LAENGE : RVECTOR[0..N-1]; HR : REAL; BEGIN FOR I := 0 TO N - 1 DO LAENGE[I] := SQR( DIAM( REST[I] ) ) * ( RVECTOR( MRA[*,I] ) * RVECTOR( MRA[*,I] ) ); FOR J := 0 TO N - 2 DO FOR K := J + 1 TO N - 1 DO IF LAENGE[J] < LAENGE[K] THEN BEGIN HR := LAENGE[K]; LAENGE[K] := LAENGE[J]; LAENGE[J] := HR; HV := RVECTOR(MRA[*,K]); MRA[*,K] := MRA[*,J]; MRA[*,J] := HV; END; QR( MRA ); (* MRA WIRD MIT Q AUS DER QR-ZERLEGUNG VON MRA ÜBERSCHRIEBEN *) END; (*--------------------------------------------------------------------------*) PROCEDURE IVP_STEP( VAR X_AKT : REAL; H : REAL; VARIANTEN: INTEGER ); (* HERZSTÜCK DES HAUPTPROGRAMMS: *) (* AUSFÜHREN EINES INTEGRATIONSSCHRITTS DES IVPS *) (* IN DIESER PROZEDUR WERDEN DIE BENÖTIGTEN DATEN *) (* FÜR DAS UNTERPROGRAMM LIVP_MPS BEREITGESTELLT *) (* DIE EIGENTLICHE INTEGRATION ERFOLGT IM UNTERPROGRAMM *) VAR I , J , K , M : INTEGER; MP_TRP : MPVECTOR[0..KMAX]; VRY : RVECTOR[0..N-1]; (*----------------*) FUNCTION RECESS( K : INTEGER ) : BOOLEAN; (* SCHÄTZT AB, OB DIE KOEFFIZIENTEN AB DEM K-TEN BETRAGSMÄSSIG FALLEN *) (* OBACHT: RECESS HAT SEITENEFFEKTE! *) (* DER AUFRUF VON CAUCHYTEST BESETZT DIE VARIABLE RMAX_LB *) VAR I , J , NU , MU : INTEGER; MINORD , MAXORD : INTEGER; RADIUS , RMAX_ALT : REAL; KOEFSUM1 , KOEFSUM2 : REAL; KOEFSUM3 , KOEFSUM : REAL; LESSONE : BOOLEAN; ABBRUCH : INTERVAL; BEGIN (* CAUCHY-SCHRANKEN AUSRECHNEN *) ERR := 0; LESSONE := FALSE; RADIUS := H / OMEGA; IF RADIUS <= 10 ** ( 300 / KMAX ) THEN BEGIN SEKTOREN := 32; ORDNUNG := 15; MINORD := 10; MAXORD := 50; IF BSPNR = 3 THEN BEGIN IF INHOM THEN (* NUR INHOMOGENITÄT MACHT PROBLEME *) CAUCHYTEST( BSPNR , ALPHA , BETA , GAMMA , X_AKT , RADIUS , ORDNUNG , TIB[-1] , SEKTOREN , SEKTOR_MAX , RMAX_MID_LB , RMAX_LB , ERR , OUTPUT ) ELSE BEGIN (* SETZEN EINIGER DUMMY-WERTE, DAMIT DIE SCHLEIFE VERLASSEN WIRD *) ERR := 0; RMAX_MID_LB := 0; RMAX_LB := 0; MAXORD := ORDNUNG; MINORD := ORDNUNG; END; END ELSE CAUCHYTEST( BSPNR , ALPHA , BETA , GAMMA , X_AKT , RADIUS , ORDNUNG , TIB[0] , SEKTOREN , SEKTOR_MAX , RMAX_MID_LB , RMAX_LB , ERR , OUTPUT ); IF ERR = 0 THEN BEGIN IF RMAX_MID_LB <= RMAX_LB / 100 THEN MAXORD := ORDNUNG ELSE MINORD := ORDNUNG; END; END; IF ( ( RADIUS > 10 ** ( 300 / KMAX ) ) OR ( ERR > 0 ) ) THEN BEGIN (* FUNKTIONSWERTE NICHT BERECHENBAR, SCHRITTWEITE VERKLEINERN *) RECESS := FALSE; WRITELN( OUTFILE ); WRITE( OUTFILE, 'Funktionswerte auf Kreisrand mit Radius ' ); WRITELN( OUTFILE, RADIUS , ' nicht berechenbar.' ); WRITELN( OUTFILE ); FLUSH( OUTFILE ); WRITELN( CONSOLE ); WRITE( CONSOLE, 'Funktionswerte auf Kreisrand mit Radius ' ); WRITELN( CONSOLE, RADIUS , ' nicht berechenbar.' ); WRITELN( CONSOLE ); FLUSH( CONSOLE ); END ELSE BEGIN (* ORDNUNG DER TAYLOR-APPROXIMATION OPTIMIEREN *) WHILE MINORD < MAXORD - 1 DO BEGIN (* ORDUNG BESTIMMEN *) IF MINORD <= MAXORD - 10 THEN ORDNUNG := ORDNUNG + 5 ELSE ORDNUNG := ( MINORD + MAXORD ) DIV 2; IF BSPNR = 3 THEN (* NUR INHOMOGENITÄT MACHT PROBLEME *) CAUCHYTEST( BSPNR , ALPHA , BETA , GAMMA , X_AKT , RADIUS , ORDNUNG , TIB[-1] , SEKTOREN , SEKTOR_MAX , RMAX_MID_LB , RMAX_LB , ERR , OUTPUT ) ELSE CAUCHYTEST( BSPNR , ALPHA , BETA , GAMMA , X_AKT , RADIUS , ORDNUNG , TIB[0] , SEKTOREN , SEKTOR_MAX , RMAX_MID_LB , RMAX_LB , ERR , OUTPUT ); IF ( RMAX_MID_LB <= RMAX_LB / 100 ) THEN MAXORD := ORDNUNG ELSE MINORD := ORDNUNG; END; ORDNUNG := MAXORD; (* SCHRANKEN B_I = VRB[I] BERECHNEN *) (* BEISPIEL 1: Y" = ALPHA e^x Y + e^(-x) - ALPHA *) IF BSPNR = 1 THEN BEGIN VRM[1] := -1; VRB[1] := 0; VRM[0] := ORDNUNG; VRB[0] := RMAX_LB * ( ( H / OMEGA ) ** N ); VRM[-1] := ORDNUNG; VRB[-1] := RMAX_LB * ( ( H / OMEGA ) ** N ); END; (* BEISPIEL 2: Y" = - ( ALPHA + BETA * COS( GAMMA*X) ) Y *) IF BSPNR = 2 THEN BEGIN VRM[1] := -1; VRB[1] := 0; VRM[0] := ORDNUNG; VRB[0] := RMAX_LB * ( ( H / OMEGA ) ** N ); END; (* BSP 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 VRM[1] := -1; VRB[1] := 0; VRM[0] := 1; VRB[0] := 0; VRM[-1] := ORDNUNG; VRB[-1] := RMAX_LB * ( ( H / OMEGA ) ** N ); 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 VRM[3] := ORDNUNG; VRB[3] := 6 *> ALPHA *> FAC( ORDNUNG ) *> RMAX_LB *> ( ( H /> OMEGA ) ** N ); (* WEITERE CAUCHY-SCHRANKEN AUSRECHNEN *) CAUCHYTEST( 41 , ALPHA , BETA , GAMMA , X_AKT , RADIUS , ORDNUNG , TIB[0] , SEKTOREN , SEKTOR_MAX , RMAX_MID_LB , RMAX_LB , ERR , OUTPUT ); VRM[0] := ORDNUNG; VRB[0] := BETA *> FAC( ORDNUNG + 1 ) *> RMAX_LB *> ( ( H /> OMEGA ) ** N ); VRM[1] := -1; VRB[1] := 0; VRM[2] := ORDNUNG; VRB[2] := 6 *> ALPHA *> ALPHA *> VRB[0] /> BETA; END; (* REZESSION DER TAYLORKOEFFIZIENTEN DER LÖSUNG DES IVPS TESTEN *) (* KOEFSUM WIRD APPROXIMIERT, UM DIE RECHNUNG ZU VEREINFACHEN *) KOEFSUM1 := 0; KOEFSUM2 := 0; KOEFSUM3 := 0; FOR I := 0 TO N - 1 DO FOR J := 0 TO VRM[I] DO KOEFSUM1 := KOEFSUM1 + ABS_SUP(TIB[I,J]) * ( RADIUS**( N-I+J ) ) * ( ( K - J + I + 0.0 ) ** I ) / ( ( K+1.0 ) ** N ); FOR I := 0 TO N - 1 DO IF VRM[I] >= N THEN KOEFSUM2 := KOEFSUM2 + VRB[I] * ( (K-VRM[I]+I+0.0) ** (I+1 ) ) / ( ( I + 1 ) * ( RADIUS ** I ) * ( (K+1.0) ** N ) ); FOR NU := K + 1 TO K + N DO KOEFSUM2 := KOEFSUM2 / NU; IF INHOM THEN BEGIN KOEFSUM3 := VRB[-1] / ( ( K + 1.0 ) ** N ); END; KOEFSUM := KOEFSUM1 + KOEFSUM2 + KOEFSUM3; WRITELN( 'KOEFSUM1 nach ',K:1,' Koeffizienten ca. ',KOEFSUM1 ); WRITELN( 'KOEFSUM2 nach ',K:1,' Koeffizienten ca. ',KOEFSUM2 ); WRITELN( 'KOEFSUM3 nach ',K:1,' Koeffizienten ca. ',KOEFSUM3 ); WRITELN( 'KOEFSUM nach ',K:1,' Koeffizienten ca. ',KOEFSUM ); WRITELN; WRITELN( OUTFILE, 'KOEFSUM1 nach ',K:1,' Koeffizienten ca. ',KOEFSUM1 ); WRITELN( OUTFILE, 'KOEFSUM2 nach ',K:1,' Koeffizienten ca. ',KOEFSUM2 ); WRITELN( OUTFILE, 'KOEFSUM3 nach ',K:1,' Koeffizienten ca. ',KOEFSUM3 ); WRITELN( OUTFILE, 'KOEFSUM nach ',K:1,' Koeffizienten ca. ',KOEFSUM ); WRITELN( OUTFILE ); WRITELN( CONSOLE, 'KOEFSUM1 nach ',K:1,' Koeffizienten ca. ',KOEFSUM1 ); WRITELN( CONSOLE, 'KOEFSUM2 nach ',K:1,' Koeffizienten ca. ',KOEFSUM2 ); WRITELN( CONSOLE, 'KOEFSUM3 nach ',K:1,' Koeffizienten ca. ',KOEFSUM3 ); WRITELN( CONSOLE, 'KOEFSUM nach ',K:1,' Koeffizienten ca. ',KOEFSUM ); WRITELN( CONSOLE ); FLUSH( CONSOLE ); IF KOEFSUM < 1 THEN LESSONE := TRUE; END; RECESS := LESSONE; END; (* RECESS *) (*-------------------------------------------------------------------------*) BEGIN (* IVP_STEP *) (* TAYLORKOEFFIZIENTEN DER KOEFFIZIENTENFUNKTIONEN AUSRECHNEN *) WRITELN( OUTFILE , 'Taylorkoeffizienten werden berechnet.' ); FLUSH( OUTFILE ); TAYKOEFF_P_I( BSPNR , ALPHA , BETA, GAMMA , TIB , X_AKT , KOEFPREC ); WRITELN( OUTFILE , 'Taylorkoeffizienten sind berechnet.' ); WRITELN; FLUSH( OUTFILE ); (* PRÜFEN, OB SUMMANDEN AB DEM P-TEN FALLEN *) (* OBACHT: FUNKTION RECESS HAT SEITENEFFEKTE! *) (* DER DORTIGE AUFRUF VON CAUCHYTEST BESETZT VARIABLE RMAX_LB *) IF RECESS( KAPPA ) OR ( SCHRITTE = 1 ) THEN BEGIN (* CAUCHY-SCHRANKEN AUSRECHNEN *) IF BSPNR = 3 THEN BEGIN IF INHOM THEN (* NUR INHOMOGENITÄT MACHT PROBLEME *) CAUCHY( BSPNR , ALPHA , BETA , GAMMA , X_AKT , H / OMEGA , ORDNUNG , TIB[-1] , SEKTOR_MAX , RMAX , ERR , OUTPUT ) ELSE RMAX := 0; END ELSE CAUCHY( BSPNR , ALPHA , BETA , GAMMA , X_AKT , H / OMEGA , ORDNUNG , TIB[0] , SEKTOR_MAX , RMAX , ERR , OUTPUT ); (* BEISPIEL 1: Y" = ALPHA e^x Y + e^(-x) - ALPHA *) IF BSPNR = 1 THEN BEGIN VRM[1] := -1; VRB[1] := 0; VRM[0] := ORDNUNG; VRB[0] := RMAX * ( ( H / OMEGA ) ** N ); VRM[-1] := ORDNUNG; VRB[-1] := RMAX * ( ( H / OMEGA ) ** N ); END; (* BEISPIEL 2: Y" = - ( ALPHA + BETA * COS( GAMMA*X) ) Y *) IF BSPNR = 2 THEN BEGIN VRM[1] := -1; VRB[1] := 0; VRM[0] := ORDNUNG; VRB[0] := RMAX * ( ( H / OMEGA ) ** N ); END; (* BSP 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 VRM[1] := -1; VRB[1] := 0; VRM[0] := 1; VRB[0] := 0; VRM[-1] := ORDNUNG; VRB[-1] := RMAX * ( ( H / OMEGA ) ** N ); 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 VRM[3] := ORDNUNG; VRB[3] := 6 *> ALPHA *> FAC( ORDNUNG ) *> RMAX *> ( ( H /> OMEGA ) ** N ); (* WEITERE CAUCHY-SCHRANKEN AUSRECHNEN *) CAUCHY( 41 , ALPHA , BETA , GAMMA , X_AKT , H / OMEGA , ORDNUNG , TIB[0] , SEKTOR_MAX , RMAX , ERR , OUTPUT ); VRM[0] := ORDNUNG; VRB[0] := BETA *> FAC( ORDNUNG + 1 ) *> RMAX *> ( ( H /> OMEGA ) ** N ); VRM[1] := -1; VRB[1] := 0; VRM[2] := ORDNUNG; VRB[2] := 6 *> ALPHA *> ALPHA *> VRB[0] /> BETA; END; IF ( SCHRITTE = 1 ) AND ( AW_REAL < 3 ) THEN BEGIN (* EINZELNEN INTEGRATIONSSCHRITT FÜR IVP DURCHFÜHREN *) LIVP( N , INHOM , VRM, VRB, TIB, KOEFPREC, MID(VIY), H, OMEGA , E_ABS , T_IND , FAK , PRODINDEX , PROD , IPROD , VIZ , OUTFILE , CONSOLE ); VIY := VIZ; WRITELN; WRITELN( 'Einschließungen an der Stelle ' , H , ':' ); WRITELN; WRITE( 'Ordnung der Polynomanteile der Taylorreihen der '); WRITELN( 'Koeffizientenfunktionen: ' , ORDNUNG : 1 ); WRITE( 'Cauchysche Koeffizientenabschätzung des Reihenrests : ' ); WRITELN( RMAX : 1 ); WRITELN; FOR I := 0 TO N - 1 DO BEGIN WRITE( 'y^(' , I : 1 , ') (' , H : 10 : 9 , ') î ' ); IWRITELN( VIY[I] ); END; END ELSE BEGIN (* IF FALSE( ( SCHRITTE = 1 ) AND ( AW_REAL < 3 ) ) *) (* INTEGRATIONSSCHRITT IM "MEHRSCHRITT-VERFAHREN" DURCHFÜHREN *) X_AKT := X_AKT + H; WRITELN( CONSOLE , 'X_AKT = ' , X_AKT ); FLUSH( CONSOLE ); WRITELN; WRITELN( 'Einschließungen an der Stelle ' , X_AKT , ':' ); WRITELN; WRITE( 'Ordnung der Polynomanteile der Taylorreihen der '); WRITELN( 'Koeffizientenfunktionen: ' , ORDNUNG : 1 ); WRITELN( 'Cauchysche Koeffizientenabschätzung des Reihenrests : ' ); WRITELN( RMAX : 1 ); WRITELN; (* VARIANTE 1: MITTELWERTVERFAHREN *) IF VARIANTEN > 15 THEN BEGIN VARIANTEN := VARIANTEN - 16; (* MITTELWERTLÖSUNG *) LIVP( N , INHOM , VRM , VRB , TIB , KOEFPREC , VR_PART1 , H , OMEGA , E_ABS , T_IND , FAK , PRODINDEX , PROD , IPROD , VI_PART1 , OUTFILE , CONSOLE ); IF VIY1 <> VINULL THEN (* INTERVALLAUSWERTUNG FÜR HOMOGENE DGL *) FOR J := 0 TO N - 1 DO BEGIN VRY := MRID[J]; LIVP( N , FALSE , VRM , VRB , TIB , KOEFPREC , VRY , H , OMEGA , E_ABS , T_IND , FAK , PRODINDEX , PROD , IPROD , VIZ , OUTFILE , CONSOLE ); MIA[*,J] := VIZ; END; VIY1 := MIA * VIY1 + VI_PART1; WRITELN( 'Mittelwertverfahren:' ); WRITELN; FOR I := 0 TO N - 1 DO BEGIN WRITE( 'y^(' , I : 1 , ') (' , X_AKT : 10 : 9 , ') î ' ); IWRITELN( VIY1[I] ); END; VR_PART1 := MID( VIY1 ); VIY1 := VIY1 - VR_PART1; WRITELN; END; (* VARIANTE 2: PARALELLEPIPEDE *) (* MITTELPUNKTSMATRIX ALS FUNDAMENTALSYSTEM *) (* INHOMOGENE DGL MIT HOMOGENEN ABD *) (* FÜR INHOMOGENE DGLN FRAGWÜRDIG *) IF VARIANTEN > 7 THEN BEGIN VARIANTEN := VARIANTEN - 8; IF INHOM THEN (* PARTIKULÄRLÖSUNG MIT HOMOGENEN ANFANGSWERTEN *) LIVP( N , TRUE , VRM , VRB , TIB , KOEFPREC , VRNULL , H , OMEGA , E_ABS , T_IND , FAK , PRODINDEX , PROD , IPROD , VI_INH , OUTFILE , CONSOLE ); (* INTERVALLAUSWERTUNG FÜR HOMOGENE DGL *) FOR J := 0 TO N - 1 DO BEGIN VRY := MRB2[*,J]; LIVP( N , FALSE , VRM , VRB , TIB , KOEFPREC , VRY , H , OMEGA , E_ABS , T_IND , FAK , PRODINDEX , PROD , IPROD , VIZ , OUTFILE , CONSOLE ); MIA[*,J] := VIZ; END; MRB2 := MID( MIA ); INV( MRB2 , MIBINV2 , K ); VIU2 := ( MIBINV2 * MIA ) * VIU2 + MIBINV2 * VI_INH; VIY2 := MRB2 * VIU2; MIWRITE( OUTFILE , 'MIA' , MIA ); MRWRITE( OUTFILE , 'MRB2' , MRB2 ); MIWRITE( OUTFILE , 'MIBINV2' , MIBINV2 ); MIWRITE( OUTFILE , 'MIBINV2 * MIA' , MIBINV2 * MIA ); VIWRITE( OUTFILE , 'VIU2' , VIU2 ); WRITELN( 'Fortsetzung mit Parallelepipeden:' ); WRITELN; FOR I := 0 TO N - 1 DO BEGIN WRITE( 'y^(' , I : 1 , ') (' , X_AKT : 10 : 9 , ') î ' ); IWRITELN( VIY2[I] ); END; WRITELN; END; (* VARIANTE 3: QR-ZERLEGUNG *) (* INHOMOGENE DGL MIT HOMOGENEN ABD *) (* FÜR INHOMOGENE DGLN FRAGWÜRDIG *) IF VARIANTEN > 3 THEN BEGIN VARIANTEN := VARIANTEN - 4; IF INHOM THEN (* PARTIKULÄRLÖSUNG MIT HOMOGENEN ANFANGSWERTEN *) LIVP( N , TRUE , VRM , VRB , TIB , KOEFPREC , VRNULL , H , OMEGA , E_ABS , T_IND , FAK , PRODINDEX , PROD , IPROD , VI_INH , OUTFILE , CONSOLE ); (* INTERVALLAUSWERTUNG FÜR HOMOGENE DGL *) FOR J := 0 TO N - 1 DO BEGIN VRY := MRB3[*,J]; LIVP( N , FALSE , VRM , VRB , TIB , KOEFPREC , VRY , H , OMEGA , E_ABS , T_IND , FAK , PRODINDEX , PROD , IPROD , VIZ , OUTFILE , CONSOLE ); MIA[*,J] := VIZ; END; MRB3 := MID( MIA ); QR_SORT( MRB3 , VIU3 ); INV( MRB3 , MIBINV3 , K ); VIU3 := ( MIBINV3 * MIA ) * VIU3 + MIBINV3 * VI_INH; VIY3 := MRB3 * VIU3; MIWRITE( OUTFILE , 'MIA' , MIA ); MRWRITE( OUTFILE , 'MRB3' , MRB3 ); MIWRITE( OUTFILE , 'MIBINV3' , MIBINV3 ); MIWRITE( OUTFILE , 'MIBINV3 * MIA' , MIBINV3 * MIA ); VIWRITE( OUTFILE , 'VIU3' , VIU3 ); WRITELN( 'Fortsetzung mit QR-Zerlegung:' ); WRITELN; FOR I := 0 TO N - 1 DO BEGIN WRITE( 'y^(' , I : 1 , ') (' , X_AKT : 10 : 9 , ') î ' ); IWRITELN( VIY3[I] ); END; WRITELN; END; (* VARIANTE 4: ERWEITERTES MITTELWERTVERFAHREN *) (* MITTELPUNKTSMATRIX ALS FUNDAMENTALSYSTEM *) IF VARIANTEN > 1 THEN BEGIN VARIANTEN := VARIANTEN - 2; (* MITTELWERTLÖSUNG *) LIVP( N , INHOM , VRM , VRB , TIB , KOEFPREC , VR_PART4 , H , OMEGA , E_ABS , T_IND , FAK , PRODINDEX , PROD , IPROD , VI_PART4 , OUTFILE , CONSOLE ); IF VIY4 <> VINULL THEN BEGIN (* HOMOGENE DGL *) FOR J := 0 TO N - 1 DO BEGIN VRY := MRB4[*,J]; LIVP( N , FALSE , VRM, VRB, TIB , KOEFPREC , VRY , H , OMEGA , E_ABS , T_IND , FAK , PRODINDEX , PROD , IPROD , VIZ , OUTFILE , CONSOLE ); MIA[*,J] := VIZ; END; END; (* PARALLELEPIPED-AUSWERTUNG *) MRB4 := MID( MIA ); INV( MRB4 , MIBINV4 , K ); VIY4 := MIA * VIU4 + VI_PART4; WRITELN; WRITE( 'Erweitertes Mittelwertverfahren mit ' ); WRITELN( 'Mittelpunktsmatrix:' ); WRITELN; FOR I := 0 TO N - 1 DO BEGIN WRITE( 'y^(' , I : 1 , ') (' , X_AKT : 10 : 9 , ') î ' ); IWRITELN( VIY4[I] ); END; VIY4 := MRB4 * ( ( MIBINV4 * MIA ) * VIU4 + MIBINV4 * VI_PART4 ); VR_PART4 := MID( VIY4 ); VIU4 := ( MIBINV4 * MIA ) * VIU4 + MIBINV4*( VI_PART4-VR_PART4 ); WRITELN; END; (* VARIANTE 5: ERWEITERTES MITTELWERTVERFAHREN MIT QR-ZERLEGUNG *) IF VARIANTEN > 0 THEN BEGIN (* MITTELWERTLÖSUNG *) LIVP( N , INHOM , VRM , VRB , TIB , KOEFPREC , VR_PART5 , H , OMEGA , E_ABS , T_IND , FAK , PRODINDEX , PROD , IPROD , VI_PART5 , OUTFILE , CONSOLE ); IF VIY5 <> VINULL THEN BEGIN (* HOMOGENE DGL *) FOR J := 0 TO N - 1 DO BEGIN VRY := MRB5[*,J]; LIVP( N, FALSE, VRM , VRB , TIB , KOEFPREC , VRY , H , OMEGA , E_ABS , T_IND , FAK , PRODINDEX , PROD , IPROD , VIZ , OUTFILE , CONSOLE ); MIA[*,J] := VIZ; END; END; VIY5 := MIA * VIU5 + VI_PART5; WRITELN; WRITELN( 'Erweitertes Mittelwertverfahren mit QR-Zerlegung:' ); WRITELN; FOR I := 0 TO N - 1 DO BEGIN WRITE( 'y^(' , I : 1 , ') (' , X_AKT : 10 : 9 , ') î ' ); IWRITELN( VIY5[I] ); END; (* QR-AUSWERTUNG *) MRB5 := MID( MIA ); QR_SORT( MRB5 , VIU5 ); INV( MRB5 , MIBINV5 , K ); VIY5 := MRB5 * ( ( MIBINV5 * MIA ) * VIU5 + MIBINV5 * VI_PART5 ); VR_PART5 := MID( VIY5 ); VIU5 := ( MIBINV5 * MIA ) * VIU5 + MIBINV5*( VI_PART5-VR_PART5 ); WRITELN; END; END; END ELSE BEGIN (* IF FALSE( RECESS( KAPPA ) OR ( SCHRITTE = 1 ) ) *) (* ANZAHL DER TEILINTERVALLE ERHÖHEN *) SUBINTERVALS := SUBINTERVALS + 1; IVP_STEP( X_AKT , H / 2 , VARIANTEN ); IVP_STEP( X_AKT , H / 2 , VARIANTEN ); END; END; (* *) (*--------------------------------------------------------------------------*) (* *) (* IVP_HAUPT *) (* *) (*--------------------------------------------------------------------------*) (* *) BEGIN (* AWP EINLESEN *) EINLESEN; (* KOPFZEILE DRUCKEN *) WRITE_HEAD; X_AKT := 0; VRNULL := 0; VINULL := 0; MRID := ID( MRID ); (* ANFANGSWERTE FÜR SPEZIELLE AUSWERTUNGEN VORBESETZEN *) MIA := MRID; MRB := MRID; MRB2 := MRID; MRB3 := MRID; MRB4 := MRID; MRB5 := MRID; VIU2 := VIY; VIU3 := VIY; VIY0 := VIY; IF AW_REAL < 3 THEN BEGIN VR_PART1 := VRY; VIY1 := 0; END ELSE BEGIN VR_PART1 := MID( VIY ); VIY1 := VIY - VR_PART1; END; VIY2 := VIY; VIY3 := VIY; VR_PART4 := VR_PART1; VIY4 := VIY1; VR_PART5 := VR_PART1; VIY5 := VIY1; VIU4 := VIY4; VIU5 := VIY5; INIT_TIMER; FAK[0] := 1; FOR I := 1 TO 30 DO FAK[I] := FAK[I-1] * I; FOR I := 0 TO N DO FOR J := 0 TO I DO BEGIN PROD[I-J,J] := FAK[I] / FAK[I-J]; IPROD[I-J,J] := PROD[I-J,J]; END; PRODINDEX := 0; SUBINTERVALS := 1; IVP_STEP( X_AKT , X_END - X_AKT , VARIANTEN ); WRITELN; WRITE( 'Integrationsbereich unterteilt in ', SUBINTERVALS : 1 ); WRITELN( ' Teilintervalle.' ); END; (* *) (*--------------------------------------------------------------------------*) (* *) (* HAUPTPROGRAM *) (* *) (*--------------------------------------------------------------------------*) (* *) BEGIN (* AWP EINLESEN *) EINLESEN_MAIN; IVP_HAUPT( N ); CPU_TIME := GET_TIMER / 1000; WRITELN; WRITELN( 'Rechenzeit:', CPU_TIME : 7 : 2, ' Sekunden.' ); END. (*==========================================================================*) (* *) (* PASCAL-XSC - PROGRAMM LIVPTAYP *) (* *) (* PROGRAMM 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) analytisch *) (* *) (* POTENZREIHENANSATZ: y(x) = \sum_0^infty a_k x^k *) (* LÖSUNG MIT VARIABLER TAYLORENTWICKLUNG UND RESTGLIEDABSCHÄTZUNG *) (* RECHNUNG MIT ADAPTIVER GENAUIGKEIT *) (* *) (* HIER : HAUPTPROGRAMM *) (* EINLESEN DES ANFANGSWERTPROBLEMS *) (* ZERLEGEN DES DEFINITIONSINTERVALLS IN GEEIGNETE ABSCHNITTE *) (* AUFRUF DER PROZEDUR LIVP IM MODUL LIVP_MPS *) (* *) (*==========================================================================*)