XProfan

Titre de la fenêtre " Anfangswert-bestimmte Differentialgleichung "+\
"1. Ordre avec PRÄDIKTOR-KORREKTOR-METHODE lösen"
'{ (T) Translated to XProfan-11 2014-08 by P.Specht, vienne
'***********************************************************
'*        SOLVING DIFFERENTIAL EQUATIONS OF ORDER 1        *
'*             (Prediction-correction method)              *
'* ------------------------------------------------------- *
' Beisp. pour une DGL 1. Ordre: Etwa cela Zerfallsgesetz  *
' d.radioactivité. suivante: 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²) à partir de x=0 to x=10  *
'*  with initial condition y(0)=1)                         *
'*                                                         *
'*  Contribution x begin: 0                                       *
'*  Contribution x end  : 10                                      *
'*  Contribution y begin: 1                                       *
'*  Contribution number of points: 10                             *
'*  Contribution 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
Fenêtre Style 24:Font 2:randomize:Fenêtre 0,0-%maxx,%maxy-40
declare tx![100],ty![100],x!,xi!,xf!,yi!,  fi&,n&,i&

si 1'  0=Manual input, 1=Example

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

d'autre

    imprimer "\n\n Contribution x begin: ";  :input xi!
    imprimer " Contribution x end  : ";      :input xf!
    Imprimer " Contribution y begin: ";      :input yi!
    imprimer " Contribution number of points: ";:input n&
    imprimer " Contribution Granularity: ";  :input fi&

endif

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

WhileLoop n&:i&=&Boucle

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

endwhile

imprimer " ----------------------------- "
waitinput
FIN
'} Fin of main program
' ici qui définition qui Differentialgleichung f:
' f: y'=( 4xy + 4x*sqrt(y))/(1+x^2)

proc f :parameters x!,y!

    declare f!

    si 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 l'ordre 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 le following steps:                 *
'*     1. calculate y2, y3, y4 using a classical Runge-   *
'*        Kutta method of l'ordre 4.                        *
'*     2. à partir de 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 le y'(n+1) = f(x(n+1),y(n+1)) by le    *
'*        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! : cas (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&=&Boucle : ni&=(i&-1)*fi&-1

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

            si k&<4

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

            d'autre

                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