XProfan

Título de la ventana " Anfangswert-cierto Differentialgleichung "+\
"1. Orden con PRÄDIKTOR-KORREKTOR-METHODE lösen"
'{ (T) Translated to XProfan-11 2014-08 by P.Pájaro carpintero, Wien
'***********************************************************
'*        SOLVING DIFFERENTIAL EQUATIONS OF ORDER 1        *
'*             (Prediction-correction method)              *
'* ------------------------------------------------------- *
' Beisp. para una DGL 1. Orden: Etwa el 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)                         *
'*                                                         *
'*  Entrada x begin: 0                                       *
'*  Entrada x end  : 10                                      *
'*  Entrada y begin: 1                                       *
'*  Entrada number of points: 10                             *
'*  Entrada 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
Ventana de Estilo 24:Font 2:randomize:Ventana 0,0-%maxx,%maxy-40
declarar tx![100],ty![100],x!,xi!,xf!,yi!,  fi&,n&,i&

if 1'  0=Manual input, 1=Example

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

más

    imprimir "\n\n Entrada x begin: ";  :input xi!
    imprimir " Entrada x end  : ";      :input xf!
    Imprimir " Entrada y begin: ";      :input yi!
    imprimir " Entrada number of points: ";:input n&
    imprimir " Entrada Granularity: ";  :input fi&

endif

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

WhileLoop n&:i&=&Loop

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

endwhile

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

proc f :parámetros x!,y!

    declarar f!

    if y!>=0

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

    endif

    volver f!

ENDPROC

proc runge_kutta :parámetros h!,x!,y!

    ' classical Runge-Kutta method of order 4
    declarar 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 :parámetros tx![],ty![],xi!,xf!,yi!,m&,fi&

    declarar ni&,h!,w!,x!,y!,z!,p![3] , k&,i&,j&
    z!=yi! : caso (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!'<<< se en Fortran90 en runge_kutta incluso como Rückparameter hecho!
                p![3-k&]=f(x!,z!)

            más

                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 bucle

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

    endwhile'i bucle

ENDPROC