XProfan

Italia

p.specht

Ist die Differentialgleichung lediglich von 1. Ordnung, determinieren also z.B. nur aktuelle Bestands-Augenblickswerte die Änderungsgeschwindigkeit eines Systems (Beispiel: Radioaktiver Zerfall), dann bietet das Prädiktor-Korrektor-Verfahren eine rasch konvergierende Lösung zur schrittweisen Berechnung der Systemverhaltens. Bestimmte Anfangsbedingungen müssen allerdings vorgegeben werden. Man spricht vom "Anfangswertproblem", das deutlich leichter lösbar ist als die Frage, welche Funktion stets bestimmte Randwerte/Randbedinungen erfüllt (Randwertproblem, oft besser mit analytischen Überlegungen lösbar - manchmal aber auch garnicht).

Nachstehend wieder eine Rückübersetzung als DEMO OHNE GEWÄHR, aus Fortran-90 bzw. ursprünglich Turbopascal. F90 bedient sich bei Subroutinen des Call by Reference, sodaß Parameter-Arithmetik sich auf das aufrufende Programm auswirkt (Stolperfalle!). Bitte nicht wundern: Für den Start-Schritt necessario das Programm den einfachen Runge-Kutta-Algorithmus.

Kompilieren Markieren Separieren

WindowTitle " Anfangswert-bestimmte Differentialgleichung "+\
"1. Ordnung mit PRÄDIKTOR-KORREKTOR-METHODE lösen"
'{ (T) Translated to XProfan-11 2014-08 by P.Specht, Wien
'***********************************************************
'*        SOLVING DIFFERENTIAL EQUATIONS OF ORDER 1        *
'*             (Prediction-correction method)              *
'* ------------------------------------------------------- *
' Beisp. per eine DGL 1. Ordnung: Etwa das Zerfallsgesetz  *
' d.Radioaktivität. Nachstehend: Willkürliche Testfunktion *
'* ------------------------------------------------------- *
'* Reference:                                              *
'*   "Analyse en Turbo Pascal versions 5.5 et 6.0 By Marc  *
'*    DUCAMP et Alain REVERCHON - Eyrolles, Paris 1991"    *
'* F90 version by J-P Moreau, Paris   (www.jpmoreau.fr)    *
'* ------------------------------------------------------- *
'* SAMPLE RUN:                                             *
'* (Integrate y'=(4xy+4x*sqrt(y))/(1+x²) from x=0 to x=10  *
'*  with initial condition y(0)=1)                         *
'*                                                         *
'*  Input x begin: 0                                       *
'*  Input x end  : 10                                      *
'*  Input y begin: 1                                       *
'*  Input number of points: 10                             *
'*  Input finesse: 10                                      *
'*                                                         *
'*         X               Y                               *
'*  -----------------------------                          *
'*     1.000000         8.999984                           *
'*     2.000000        80.999881                           *
'*     3.000000       360.999496                           *
'*     4.000000      1088.998512                           *
'*     5.000000      2600.996484                           *
'*     6.000000      5328.992838                           *
'*     7.000000      9800.986875                           *
'*     8.000000     16640.977767                           *
'*     9.000000     26568.964559                           *
'*    10.000000     40400.946171                           *
'*  -----------------------------                          *
'*                                                         *
'}**********************************************************
'{ PROGRAM EQUDIFPC
WindowStyle 24:Font 2:randomize:Window 0,0-%maxx,%maxy-40
declare tx![100],ty![100],x!,xi!,xf!,yi!,  fi&,n&,i&

if 1'  0=Manual input, 1=Example

    xi!=0  :print "\n\n X_begin: ";xi!
    xf!=10 :print " X_end  : ";xf!
    yi!=1  :Print " Y_begin: ";yi!
    n& =10 :print " Number of points inbetween: ";n&
    fi&=10 :print " Granularity between points: ";fi&

else

    print "\n\n Input x begin: ";  :input xi!
    print " Input x end  : ";      :input xf!
    Print " Input y begin: ";      :input yi!
    print " Input number of points: ";:input n&
    print " Input Granularity: ";  :input fi&

endif

equadiff_pc tx![],ty![],xi!,xf!,yi!,n&,fi&
print' Results:
print "      X                Y      "
print " -----------------------------"

WhileLoop n&:i&=&Loop

    print format$("%g",tx![i&]),format$("%g",ty![i&])

endwhile

print " ----------------------------- "
waitinput
END
'} End of main program
' Hier die Definition der Differentialgleichung f:
' f: y'=( 4xy + 4x*sqrt(y))/(1+x^2)

proc f :parameters x!,y!

    declare f!

    if y!>=0

        f! = ( 4*x!*y!+4*x!*sqrt(y!) ) / ( 1+sqr(x!) )

    endif

    return f!

endproc

proc runge_kutta :parameters h!,x!,y!

    ' classical Runge-Kutta method of order 4
    declare a!,b!,c!,d!,f!
    a!=h!*f(x!,y!)
    b!=h!*f(x!+h!/2, y!+a!/2)
    c!=h!*f(x!+h!/2, y!+b!/2)
    d!=h!*f(x!+h!,y!+c!)
    y!=y!+(a!+2*b!+2*c!+d!)/6
    z!=y!

endproc

'**********************************************************
'{             Prediction-correction method               *
'* ------------------------------------------------------ *
'* INPUTS:   xi      begin x value                        *
'*	         xf      end x value                          *
'*           y1      begin y value (at x=xi)              *
'*           m       number of points to calculate        *
'*           fi      finesse (number of intermediate      *
'*                   points (for example 10)              *
'* OUTPUTS:                                               *
'*           tx      table of x values (m values)         *
'*           ty      table of y values (m values)         *
'*                                                        *
'* DESCRIPTION:                                           *
'* The algorithm has the following steps:                 *
'*     1. calculate y2, y3, y4 using a classical Runge-   *
'*        Kutta method of order 4.                        *
'*     2. from point (x4, y4), first estimate y(n+1) by   *
'*        PREDICTION formula:                             *
'*  y(n+1)=y(n)+h/24(55y'(n)-59y'(n-1)+37y'(n-2)-9y'(n-3) *
'*                                                        *
'*  then continue with CORREKTION formula:                *
'*  y(n+1)=y(n) + h/24(9y'(n+1)+19y'(n)-5y'(n-1)+y'(n-2)  *
'*                                                        *
'*  substituting the y'(n+1) = f(x(n+1),y(n+1)) by the    *
'*        estimated value of y(n+1), until convergence    *
'*        is obtained.                                    *
'}*********************************************************

Proc equadiff_pc :parameters tx![],ty![],xi!,xf!,yi!,m&,fi&

    declare ni&,h!,w!,x!,y!,z!,p![3] , k&,i&,j&
    z!=yi! : case (m&>100) or (fi&<1) : Goto "retour"
    h!=(xf!-xi!)/fi&/m&
    p![3]=f(xi!,yi!) : tx![0]=xi! : ty![0]=yi!
    k&=0

    whileloop m&:i&=&Loop : ni&=(i&-1)*fi&-1

        whileloop fi&:j&=&Loop : x!=xi!+h!*(ni&+j&) : inc k&

            if k&<4

                runge_kutta h!,x!,z!
                x!=x!+h!'<<< wird in Fortran90 in runge_kutta selbst als Rückparameter erledigt!
                p![3-k&]=f(x!,z!)

            else

                x!=x!+h!
                w!=z! +h!/24*(55*p![0]-59*p![1] +37*p![2]  -9*p![3])
                ' y(n+1)=y(n)+h/24*(55*y'(n)-59y'(n-1)+37*y'(n-2)-9*y'(n-3)

                Repeat

                    y!=w!
                    w!=z!+h!/24*(9*f(x!,y!)+19*p![0]-5*p![1] + p![2])
                    '   y(n+1)=y(n)+h/24(9*y'(n+1) +19*y'(n)-5*y'(n-1)+y'(n-2)

                Until abs(y!-w!)<val("1e-10")

                z!=w! : p![3]=p![2] : p![2]=p![1] : p![1]=p![0]
                p![0]=f(x!,z!)

            endif

        endwhile'j loop

        tx![i&]=x!
        ty![i&]=z!

    endwhile'i loop

endproc