XProfan

English

p.specht

is The Differentialgleichung solely of 1. Order, determinieren means z.B. only actually Bestands-Augenblickswerte The Änderungsgeschwindigkeit one Systems (example: Radioaktiver Zerfall), then bid the Prädiktor-Korrektor-take action a rasch konvergierende Solution to schrittweisen Berechnung the Systemverhaltens. defined Anfangsbedingungen must though pretended go. one speaks of "Anfangswertproblem", the explicit leichter lösbar is as The question, which function always defined Randwerte/Randbedinungen erfüllt (Randwertproblem, often rather with analytischen Überlegungen lösbar - sometimes but too not at all).

hereinafter again a Rückübersetzung as DEMO OHNE GEWÄHR, from Fortran-90 or. original Turbopascal. F90 bedient itself with Subroutinen the Call by Reference, sodaß Parameter-Arithmetik itself the aufrufende Program auswirkt (Stolperfalle!). Please not wonder: for the Start-step needed the program whom einfachen Runge-Kutta-Algorithmus.

Compile Mark Separation

Window Title " Anfangswert-defined Differentialgleichung "+\
"1. Order with PRÄDIKTOR-KORREKTOR-METHODE lösen"
'{ (T) Translated to XProfan-11 2014-08 by P.woodpecker, Wien
'***********************************************************
'*        SOLVING DIFFERENTIAL EQUATIONS OF ORDER 1        *
'*             (Prediction-correction method)              *
'* ------------------------------------------------------- *
' Beisp. for a DGL 1. Order: about the Zerfallsgesetz  *
' d.radioactivity. hereinafter: 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
Window Style 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 stature$("%g",tx![i&]),stature$("%g",ty![i&])

endwhile

print " ----------------------------- "
waitinput
END
'} End of main program
' here The Definition the 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), ridge 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!'<<< becomes in Fortran90 in runge_kutta self as Rückparameter Done!
                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