Monday, September 18, 2017

Adams-Moulton Prediction-Correction Method

'********************************************************************
'*   Solve Y' = F(X,Y) with initial conditions using the Adams-     *
'*   Moulton Prediction-Correction Method                           *
'* ---------------------------------------------------------------- *
'* SAMPLE RUN                                                       *
'* (Integrate Y' = -Y + X/(1+X)^2 from X=0 to X=1 with initial      *
'*  condition Y(0) = 1 )                                            *
'*                                                                  *
'*     X          Y         Y True   |Y-Y True|                     *
'* ---------------------------------------------                    *
'*  0.000000   1.000000   1.000000    0.000000                      *
'*  0.050000   0.952381   0.952381    0.000000                      *
'*  0.100000   0.909091   0.909091    0.000000                      *
'*  0.150000   0.869569   0.869565    0.000004   1                  *
'*  0.200000   0.833340   0.833333    0.000006   1                  *
'*  0.250000   0.800008   0.800000    0.000009   1                  *
'*  0.300000   0.769241   0.769231    0.000010   1                  *
'*  0.350000   0.740752   0.740741    0.000011   1                  *
'*  0.400000   0.714298   0.714286    0.000012   1                  *
'*  0.450000   0.689668   0.689655    0.000012   1                  *
'*  0.500000   0.666679   0.666667    0.000013   1                  *
'*  0.550000   0.645174   0.645161    0.000013   1                  *
'*  0.600000   0.625013   0.625000    0.000013   1                  *
'*  0.650000   0.606073   0.606061    0.000013   1                  *
'*  0.700000   0.588248   0.588235    0.000013   1                  *
'*  0.750000   0.571441   0.571429    0.000013   1                  *
'*  0.800000   0.555568   0.555556    0.000012   1                  *
'*  0.850000   0.540553   0.540541    0.000012   1                  *
'*  0.900000   0.526327   0.526316    0.000012   1                  *
'*  0.950000   0.512832   0.512821    0.000011   1                  *
'*  1.000000   0.500011   0.500000    0.000011   1                  *
'* ---------------------------------------------------------------- *
'* REFERENCE: "Méthode de calcul numérique- Tome 2 - Programmes en  *
'*             Basic et en Pascal By Claude Nowakowski, Edition du  *
'*             P.S.I., 1984".                                       *
'*                                                                  *
'*                      Quick Basic Release By J-P Moreau, Paris.   *
'*                                  (www.jpmoreau.fr)               *
'********************************************************************
'See explanation file adambash.txt
'---------------------------------
DefInt I-N
DefDbl A-H, O-Z

Option Base 0        'index begins from zero

Dim X(3), Y(3)

H = 0.05             'integration step
X(0) = 0#: Y(0) = 1# 'Initial conditions
EC = 0.000001        'Precision

F$ = " ##.######  ##.######  ##.######   ##.######  ##"

Cls
Print "     X          Y         Y True   |Y-Y True| "
Print " ---------------------------------------------"
Print USING; F$; X(0); Y(0); Y(0); X(0)

'Start with Runge-Kutta
For K = 0 To 1
  XX = X(K): YY = Y(K): GoSub 1000: C1 = C
  XX = X(K) + H / 2#: YY = Y(K) + H / 2# * C1: GoSub 1000: C2 = C
  YY = Y(K) + H / 2# * C2: GoSub 1000: C3 = C
  X(K + 1) = X(K) + H
  XX = X(K + 1): YY = Y(K) + H * C3: GoSub 1000: C4 = C
  Y(K + 1) = Y(K) + H * (C1 + 2 * C2 + 2 * C3 + C4) / 6#
  XX = X(K + 1): GoSub 2000: ER = Abs(F - Y(K + 1))

  Print USING; F$; X(K + 1); Y(K + 1); F; ER

Next K

100 K = 2
XX = X(K): YY = Y(K): GoSub 1000: C1 = C
XX = X(K - 1): YY = Y(K - 1): GoSub 1000: C2 = C
XX = X(K - 2): YY = Y(K - 2): GoSub 1000: C3 = C
X(K + 1) = X(K) + H
YP = Y(K) + H / 12# * (23 * C1 - 16 * C2 + 5 * C3)

L = 0
200 XX = X(K + 1): YY = YP: GoSub 1000: C1 = C
XX = X(K): YY = Y(K): GoSub 1000: C2 = C
XX = X(K - 1): YY = Y(K - 1): GoSub 1000: C3 = C
YC = Y(K) + H / 12# * (5 * C1 + 8 * C2 - C3)

'PRINT YC

If Abs(YP - YC) > EC Then
  YP = YC: L = L + 1: GoTo 200
End If

Y(K + 1) = YC: XX = X(K + 1): GoSub 2000: ER = Abs(F - Y(K + 1))

Print USING; F$; X(K + 1); Y(K + 1); F; ER; L


For K = 0 To 2
  X(K) = X(K + 1): Y(K) = Y(K + 1)
Next K

If X(2) < 1# Then
  GoTo 100
End If

End

1000 C = -YY + XX / ((1# + XX) ^ 2)
Return

2000 F = 1# / (1# + XX)
Return

Monday, September 11, 2017

Adams-Bashforth Method

'********************************************************************
'*    Solve Y' = F(X,Y) with Initial Condition Y(X0)=Y0 using the   *
'*    Adams-Bashforth Method                                        *
'* ---------------------------------------------------------------- *
'* REFERENCE: "Méthode de calcul numérique- Tome 2 - Programmes en  *
'*             Basic et en Pascal By Claude Nowakowski, Edition du  *
'*             P.S.I., 1984" [4].                                   *
'* ---------------------------------------------------------------- *
'* SAMPLE RUN:                                                      *
'* (Solve Y' = -Y + X/((1+X)*(1+X))  with Y(0)=1 for X0=0 up to     *
'*  X1=1.0, exact solution is Y = 1 / (1+X).                        *
'*                                                                  *
'*      X           Y        Y exact      Error                     *
'* -----------------------------------------------                  *
'*   0.000000    1.000000    1.000000   0                           *
'*   0.000000    1.000000    1.000000   0.000000                    *
'*   0.050000    0.952381    0.952381   0.000000                    *
'*   0.100000    0.909091    0.909091   0.000000                    *
'*   0.150000    0.869525    0.869565   0.000040                    *
'*   0.200000    0.833265    0.833333   0.000068                    *
'*   0.250000    0.799910    0.800000   0.000090                    *
'*   0.300000    0.769125    0.769231   0.000106                    *
'*   0.350000    0.740623    0.740741   0.000117                    *
'*   0.400000    0.714160    0.714286   0.000125                    *
'*   0.450000    0.689525    0.689655   0.000131                    *
'*   0.500000    0.666533    0.666667   0.000134                    *
'*   0.550000    0.645026    0.645161   0.000135                    *
'*   0.600000    0.624865    0.625000   0.000135                    *
'*   0.650000    0.605926    0.606061   0.000134                    *
'*   0.700000    0.588103    0.588235   0.000133                    *
'*   0.750000    0.571298    0.571429   0.000131                    *
'*   0.800000    0.555428    0.555556   0.000128                    *
'*   0.850000    0.540416    0.540541   0.000125                    *
'*   0.900000    0.526194    0.526316   0.000121                    *
'*   0.950000    0.512703    0.512821   0.000118                    *
'*   1.000000    0.499886    0.500000   0.000114                    *
'* -----------------------------------------------                  *
'*                                                                  *
'*                              Basic Release By J-P Moreau, Paris. *
'*                                       (www.jpmoreau.fr)          *
'********************************************************************
DefDbl A-H, O-Z
DefInt I-N

Option Base 0

  H = 0.05 'integration step

  Dim B(4), X(4), Y(4)

  'Initial conditions
  X(0) = 0#      'starting X
  X1 = 1#        'ending X
  Y(0) = 1#      'initial Y
  
  'write header
  Cls
  Print
  Print "       X           Y        Y exact      Error    "
  Print "  ------------------------------------------------"
  'write initial line
  F$ = "#####.######"
  Print USING; F$; X(0);
  Print USING; F$; Y(0);
  Print USING; F$; Y(0);
  Print "    0"
  'use Runge-Kutta to start
  For K = 0 To 1
    GoSub 1000 'call RK4
    XX = X(K + 1): GoSub 600: ER = Abs(FX - Y(K + 1))
    Print USING; F$; X(K + 1);
    Print USING; F$; Y(K + 1);
    Print USING; F$; FX;
    Print USING; F$; ER
  Next K
  'main integration loop
  While X(2) < X1
    For I = 1 To 3
      XX = X(3 - I): YY = Y(3 - I): GoSub 500
      B(I) = F
    Next I
    X(3) = X(2) + H
    Y(3) = Y(2) + H * (23# * B(1) - 16# * B(2) + 5# * B(3)) / 12#
    XX = X(3): GoSub 600: ER = Abs(Y(3) - FX)
    Print USING; F$; X(3);
    Print USING; F$; Y(3);
    Print USING; F$; FX;
    Print USING; F$; ER
    For K = 0 To 2
      X(K) = X(K + 1): Y(K) = Y(K + 1)
    Next K
  Wend
  Print "  ------------------------------------------------"

End

'User defined function Y'=F(XX,YY)
500  F = -YY + XX / ((1# + XX) * (1# + XX))
     Return

'Exact solution Y=FX(XX)
600  FX = 1# / (1# + XX)
     Return

'Runge-Kutta method to calculate first points only
1000 XX = X(K): YY = Y(K): GoSub 500: C1 = F
     XX = X(K) + H / 2#: YY = Y(K) + H / 2# * C1: GoSub 500: C2 = F
     XX = X(K) + H / 2#: YY = Y(K) + H / 2# * C2: GoSub 500: C3 = F
     XX = X(K) + H: YY = Y(K) + H * C3: GoSub 500: C4 = F
     X(K + 1) = X(K) + H
     Y(K + 1) = Y(K) + H * (C1 + 2 * C2 + 2 * C3 + C4) / 6#
     Return

'end of file adambash.bas



 EXPLANATION FILE OF PROGRAM ADAMBASH
====================================


  Ordinary Differential Equations Y' = F(x,y)
  -------------------------------------------


  Linked Steps Method
  -------------------

    We have seen that Euler, Runge-Kutta... methods can be put under the general
  form:

                y    = y  + h phi(x , y , h)
                 n+1    n          n   n

  and each point y    of the solution is only determined from previous point, y
                  n+1                                                          n

  Calculations are made for y1, y2,...,yn-1, yn, yn+1... These methods are cal-
  led "with separate steps": each step is independant from the previous one. To
  obtain a given accuracy, each step has to be divided into intermediate steps.

  Let us take an example:

  For the Runge-Kutta method of order 4, we use the formulas:

                k1 = h f(xn, yn)

                k2 = h f(xn+h/2, yn+k1/2)

                k3 = h f(xn+h/2, yn+k2/2)

                k4 = h f(xn+h, yn+k3)

                y    = y  + (1/6)(k1+2k2+2k3+k4)
                 n+1    n

  Another method consists in using the previous calculations to improve speed;
  the y   point is evaluated from y , y    and y    points.
       n+1                         n   n-1      n-2

  The general formula is:

                y    = a y + a y   + ... + a y 
                 n+1    0 n   1 n-1         k n-k

                       + h (b  f   + b f + b f    + ... + b f
                             -1 n+1   0 n   1 n-1          k n-k

  This can be written as:

                        k              k
                y    = Sum a y    + h Sum  b  f
                 n+1   j=0  j n-j     j=-1  j  n-j


    These methods are called "with linked steps" and f(x,y) evaluations are only
  made at points x0, x1, x2,..., xn.

  If b   = 0, the process is explicit; y    is directly obtained by applying the
      -1                                n+1
  formula.

  If b   <> 0, the process is implicit: we must solve an equation of the form
      -1

  y    = phi(y   ) to obtain y   .
   n+1        n+1             n+1

    For the three first points, y1, y2, y3, we have no previous points to calcu-
  late them: these methods cannot start by themselves like "with separate steps"
  methods.

    The truncation error Tn is estimated by using the Taylor formula for y(x   )
                                                                            n+j
                                 (jh)²                (jh)^p  (p)
    y(x   ) = y(x ) + jhy'(x ) + ----- y"(x ) + ... + ------ y   (x )
       n+j       n          n      2!      n            p!         n

                    + h^p eps(h)

    and also for f(x   , y(x   )):
                    n+j     n+j


    f(x   , y(x   )) = y'(x +jh)
       n+j     n+j         n                         p-1
                                 h               (jh)      (p)
                     = y'(x ) + j y"(x ) + ... + -------- y   (x ) + h^p eps(h)
                           n          n           (p-1)!        n

    Example: let us estimate the truncation error for the implicit formula:

                     y    = y  + (h/2) (f    + f )
                      n+1    n           n+1    n

                     T    = y(x   ) - y
                      n+1      n+1     n+1
                                                                    3
    y(x   ) = y(x ) + hy'(x ) + (h^2/2) y"(x ) + (h^3/3) y"'(x ) + h  eps(h)
       n+1       n         n                n                 n

                                                    2               3
    (h/2) f    = (h/2) y'(x + h) = (h/2) y'(x ) + (h /2) y"(x ) + (h /4)
           n+1             n                 n               n
                              3
                   y"'(x ) + h  eps(h)
                        n

    (h/2) f  = (h/2) y'(x )
           n             n
                                                                  
    Hence   y(x   ) - y    = y(x   ) - y  + (h/2) (f    + f ) = 
               n+1     n+1      n+1     n           n+1    n

                               3                 3
                           = (h /12) y"'(x )  + h  eps(h)
                                          n

    Here the method is of order two.


    The methods with separate steps can be integrated by the formula

                                    xn+2h
                y(x   ) - y(x   ) =  Sum  f(t,y(t)) dt
                   n+1       n-M     xn

    and by applying the Simpson's formula

                 b             b - a            a + h
                Sum (f(t) dt = ----- {f(x) + 4f(-----) + f(b)}
                 a               6                2

    Hence       y    - y  = (h/3)[f + 4f   + f   ]
                 n+2    n          n    n+1   n+2

    Here we have an implicit process.


    In a more general way, knowing y , y   , y   , we can calculate f , f   ,
                                    n   n-1   n-2                    n   n-1

  f    and approximate y' = f(x,y) by an interpolation polynomial at points  
   n-2

  x , x   , x   :
   n   n-1   n-2                      x
                                       n+1
                        y    = y    + Sum  P(x) dx 
                         n+1    n-M   x
                                       n-M

   where y    is an approximation of y(x   ) and y    is an approximation of
          n-M                           n-M       n+1

   y(x   ).
      n+1

  In the case of an implicit method, y(x   ) can be approximated by the formula:
                                        n+1
                                     x
                         p            n+1
                        y   = y    + Sum  P(x) dx
                         n+1   n-M   x
                                      n-M

                          p             p
  This allows evaluating f   = f(x   , y   ) and we can resume the interpolation
                          n+1     n+1   n+1

  step with a new polynomial P*(x). y    is calculated by the correction formula
                                     n+1

                                     x
                         c            n+1
                        y   = y    + Sum  P*(x) dx
                         n+1   n-M   x
                                      n-M

  As the points are equally spaced and indices are decreasing, xn, xn-1, xn-2...
  we can use the Newton formula with back diferences:

             Div (f ) = f - f   
                   n     n   n-1

                2
             Div (f ) = f - 2f    +f
                   n     n    n-1   n-2

             --------------------------
                k           |k|      |k|               k
             Div (f ) = f - | | f    | | f   +...+ (-1) f
                   n     n  |1|  n-1 |2|  n-2            n-k

  Hence              x                             2
         p            n+1      Div(fn)          Div (fn)
        y   = y    + Sum  [f + ------- (x-x ) + -------- (x-x )(x-x   ) + ...
         n+1   n-M   x      n    h         n      2h^2       n     n-1
                      n-M

                        k
                     Div (fn)
                   + -------- (x-x )(x-x   )...(x-x     )] dx
                      k! h^k      n     n-1        n-k+1

        Let us put u = (x-x )/h, then  du = dx/dh and
                           n                            2
                    p             1                  Div (fn)
                   y    = y    + Sum [f + Div(f )u + -------- u(u+1) + ...
                    n+1    n-M   -M    n       n        2

                        k
                     Div (fn)
                   + -------- u(u+1)(u+2)...(u+k-1) h du
                        k!

       After integration:

        p               _     _           _     2             _    k
       y    = y    + h [P f + P Div(f ) + P  Div (f ) + ... + P Div (f )]
        n+1    n-M       0 n   1     n     2       n           k      n

              _    1   1
       where  P  = -- Sum u(u+1)(u+2)...(u+j-1) du 
               j   j! -M

                j           |j|      |j|               j
       and   Div (f ) = f - | | f    | | f   +...+ (-1) f
                   n     n  |1|  n-1 |2|  n-2            n-j

                            |j|      j!
                    with    | | = ---------   (Newton's coefficients) 
                            |m|   m! (j-m)!

       Fianally   p
                 y    = y
                  n+1    n-M + h [P f  + P f    + ... + P f   ]
                                   0 n    1 n-1          k n-k

       Example: M=3, k=2

                     x                              2
        p             n+1       Div(fn)          Div (fn)
       y    = y    + Sum  [f  + ------- (x-x ) + -------- (x-x )(x-x   )] dx
        n+1    n-3   x      n     n         n      2h²        n     n-1
                      n-3

       Let us put u = (x-xn)/h, then
                                             2
        p             1                   Div (fn)
       y    = y    + Sum [f + Div(f ) u + -------- u (u+1)] h du
        n+1    n-3   -3    n       n         2

       After integration:

        p                                        2
       y    = y    + h [4f - 4 Div(f ) + (8/3)Div (f )
        n+1    n-3        n         n               n

  This leads to the Milne's formula (explicit process):

        p
       y    = y    + (4h/3) [2f - f   + 2f   ]
        n+1    n-3             n   n-1    n-2

  We do in a similar way for the implicit process. For example, with M=1 and
  k=2, the Milne's corrector formula is:

        c
       y    = y    + (h/3) [f    + 4f  + f   ]
        n+1    n-1           n+1     n    n-1

  This corresponds to the numerical integration of

       x
        n+1
       Sum  f(x) dx   by the Simpson's method. 
       x
        n-1

  For the Adams-Moulton's implicit formulas, we have:

              p          h
       k=1   y    = y  + - [3f - f   ]
              n+1    n   2    n   n-1

              c          h
             y    = y  + - [f   + f ]
              n+1    n   2   n+1   n

              p          h
       k=2   y    = y  + -- [23f - 16f   + 5f   ]
              n+1    n   12     n     n-1    n-2

              c          h
             y    = y  + -- [5f   + 8f  - f   ]
              n+1    n   12    n+1    n    n-1

              p          h
       k=3   y    = y  + -- [55f - 59f   + 17f    - 9f   ]
              n+1    n   24     n     n-1     n-2     n-3

              c          h
             y    = y  + -- [9f   + 19f  - 5f    +f   ]
              n+1    n   24    n+1     n     n-1   n-2

              p           h
       k=4   y    = y  + --- [1901f - 2984f   + 2616f    - 1274f   + 251f   ]
              n+1    n   720       n       n-1       n-2        n-3      n-4

              c           h
             y    = y  + --- [251f   + 646f  - 264f    +106f   - 19f   ]
              n+1    n   720      n+1      n       n-1      n-2     n-3


    In these formulas, the main part of what is left aside in the integration
  corresponds to the truncation error.
                                                 3
                  c                             h   "'
    Example: for y    = y  + (h/2) (f   + f ) - -- y  (ksi)
                  n+1    n           n+1   n    12

             with x  <= ksi <= x
                   n            n+1

    Another way to obtain an explicit or implicit process is using a recursive
  formula:

            y    = a y  + a y    + ... + a y
             n+1    0 n    1 n-1          k n-k

                                       + h [b  f    + b f  + ... + b y   ]
                                             -1 n+1    0 n          k n-k

    For y, the polynomial of greatest degree is:

           y(x) = 1,  y'(x) = 0
           y(x) = x,  y'(x) = 1
           y(x) = x², y'(x) = 2x
           ---------------------
                   m            m-1
           y(x) = x , y'(x) = mx

    So
           1 = a0 + a1 + ... + ak
           x = (x-h)a  + (x-2h)a  + ... + ha   + b  +b + b + ... +b
                     0          1           k-1   -1  0   1        k
           
           x² = (x-h)²a  + (x-2h)²a  + ... + ha
                       0           1           k-1

                       + 2 [xb  +(x-h)b + ... + hb
                              -1       0          k-1
           --------------------------------------------------
            l        l           l
           x  = (x-h) a  + (x-2h) a  + ... + ha
                       0           1           k-1

                             l-1          l-1
                       + l [x   b   +(x-h)   b + ... + hb
                                 -1           0          k-1

  Example: Nystroem's explicit formula:

                                        3
                       y    = y    + h Sum b  f
                        n+1    n-1     j=0  j  n-j

                            3         2
                       y = x , y' = 3x

                      3                3
         ==>  y    = x ,  y    = (x-2h)
               n+1         n-1
                                        2                2                2
                       f  = y'  = 3(x-h) , f    = 3(x-2h) , f    = 3(x-3h)
                        n     n             n-1              n-2

               3         2              2           2           2
         ==>  x  = (x-2h)  + 3h [b (x-h)  + b (x-2h)  + b (x-3h)
                                  0          1           2

  By developing and identifying:

                   b0 + b1 + b2 = 2

                   b0 + 2b1 +3b2 = 2

                   b0 + 4b1 + 9b2 = 8/3

     We find:      b0 = 7, b1 = -2, b2 = 1.


  So we have the Adams's formulas:

     Explicit, order 2:  y    = y  + h/2 (3f  - f   )
                          n+1    n          n    n-1

               order 3:  y    = y  + h/12 (23f  - 16f    + 5f   )
                          n+1    n            n      n-1     n-2

     Implicit, order 2:  y    = y  + h/2 (f   + f )
                          n+1    n         n+1   n

  and Nystroem's formulas:

     Explicit, order 2:  y    = y   + 2h f
                          n+1    n-1      n

               order 3:  y    = y   + h/3 (7f  - 2f    + f   )
                          n+1    n+1         n     n-1    n-2

     Implicit, order 3:  y    = y   + h/12 (f   + 4f -f   )
                          n+1    n-1         n+1    n  n-1

  From [BIBLI 04].
------------------------------------------------
End of file adambash.txt



Thursday, September 7, 2017

Euler-Romberg Method

'********************************************************************
'*    Solve Y' = F(X,Y) with Initial Condition Y(X0)=Y0 using the   *
'*    Euler-Romberg Method                                          *
'* ---------------------------------------------------------------- *
'* SAMPLE RUN:                                                      *
'* (Solve Y' = X*X*Y with Y(0)=1 for X0=0 up to Xn=1.1, exact       *
'*  solution is Y = Exp(X^3/3) ).                                   *
'*                                                                  *
'*    X         Y            Y       Error       Number of          *
'*          estimated      exact                subdivisions        *
'* ----------------------------------------------------------       *
'*   0.1     1.000333    1.000333  0.00000000       4               *
'*   0.2     1.002670    1.002670  0.00000001       4               *
'*   0.3     1.009041    1.009041  0.00000006       4               *
'*   0.4     1.021562    1.021562  0.00000014       4               *
'*   0.5     1.042547    1.042547  0.00000027       4               *
'*   0.6     1.074654    1.074655  0.00000086       4               *
'*   0.7     1.121125    1.121126  0.00000107       4               *
'*   0.8     1.186094    1.186095  0.00000126       4               *
'*   0.9     1.275067    1.275069  0.00000133       4               *
'*   1.0     1.395611    1.395612  0.00000114       4               *
'*   1.1     1.558410    1.558412  0.00000047       4               *
'* ----------------------------------------------------------       *
'*                                                                  *
'* Ref.: "Methodes de calcul numerique By Claude Nowakowski, Tome 2 *
'*        PSI Editions, France, 1981" [BIBLI 04].                   *
'*                                                                  *
'********************************************************************
'Program EulerRomberg


DefInt I-N
DefDbl A-H, O-Z

Option Base 0

  NMAX = 100
  H = 0.1          'initial integration step
  ER = 0.000001    'desired precision
  LA = 10          'maximum number of subdivisions
  NC = 10          'number of calculations NC = (Xn+1 - X1)/H

  Dim X(NMAX), Y(NMAX), T(20)

  'Initial conditions
  X(0) = 0#: Y(0) = 1#
  'write header
  Cls
  Print
  Print "   X         Y            Y       Error      Number of   "
  Print "         estimated      exact               subdivisions "
  Print "---------------------------------------------------------"
  'main integration loop
  For N = 0 To NC
    XC = X(N): YC = Y(N)
    XX = XC: YY = YC: GoSub 1000
    T(1) = Y(N) + H * F
    L = 1: LM = 2: ET = 1#
    While L < LA And ET >= ER
      XC = X(N): YC = Y(N)
      For J = 1 To LM
        XC = XC + H / LM
        XX = XC: YY = YC: GoSub 1000
        YC = YC + H / LM * F
      Next J
      T(L + 1) = YC: M = 1: K = L: MM = 2: ET = 1#
      If K > 1 Then
        While ET >= ER And K > 1
          T(K) = (MM * T(K + 1) - T(K)) / (MM - 1)
          ET = Abs(T(K) - T(K - 1))
          M = M + 1: K = K - 1: MM = MM * 2
        Wend
      End If
      If K = 1 Then
        L = L + 1: LM = LM * 2
      End If
    Wend
    X(N + 1) = X(N) + H: Y(N + 1) = T(K)
    XX = X(N + 1): GoSub 2000: YEX = FX
    EF = Abs(YEX - Y(N + 1))
    Print USING; "###.#"; X(N + 1);
    Print USING; "#####.######"; Y(N + 1);
    Print USING; "#####.######"; YEX;
    Print USING; "#####.########"; EF;
    Print "     "; L
  Next N
  Print "---------------------------------------------------------"

End 'of main program


'Y' = F(X,Y)
1000 'Function F(XX,YY)
  F = XX * XX * YY
Return

'Exact solution FX(XX)
2000 'Function FX(XX)
  FX = Exp(XX ^ 3 / 3)
Return

'end of file eulromb.bas

Wednesday, September 6, 2017

Explanation File for Euler-Romberg Method


                Explanation File for Euler-Romberg Method
                =========================================


          We want to solve the Cauchy's problem:


            | y' = f(x,y)
            |
            | with y(x0) = y0

        using the Euler-Romberg method.

          For that, we calculate approximate values for y1, y2,...yn,
        the exact solution is y(x1), y(x2),...y(xn), using an iterative
        algorithm based on:

                  * repeatedly apply the Euler method in the same
                    interval with an integration step halved after
                    each iteration.

                  * linearly extrapolate to obtain a better
                    approximation.

        This procedure allows getting the same precision all along the
        computation.

        Let be (xn, yn) and a starting step h0, the approximation yn+1
        is obtained by building up a table in the same way than in the
        integration Romberg method.
                                                                  0
        We start with estimating an initial value of yn+1, noted y0 by
        Euler method, using step h0:

               0
              y0 = yn + h0 f(xn,yn)

        then the new step is h0/2 and we do the same computations for
        all the points of the sub-interval:

                          0
                 h0      y0

                          1       0
           L=1   h0/2    y0      y1

                     2    2       1       0
           L=2   h0/2    y0      y1      y2

           ----------------------------------
                     k    k       k-1          1       0
           L=k   h0/2    y0      y       ...  y       y                      
                                  1            k-1     k
                           0
        The first element y0 of the column 0 is an approximation of
        y(x   ) with the step h0=x    - x . The other elements are the
           n+1                    n-1    n

        successive approximations of y(x   ) given by the Euler's formula
                                        n+1
        for h0/2, h0/4... i.e.:

                  L=1          h1=h0/2

                  x0=xn        y0=yn        f0=f(x0,y0)          

                  x1=xn+h1     y1=y0+h1f0   f1=f(x1,y1)

                                              1
                  x2=xn+2h1    y2=y1+h1f1 -> Y0 = Y2
                                             -------
                      k
        The elements Y  of the column m (m=1,2...) are obtained by using
                      m
        the linear extrapolation formula:

                            k              k+1
                    k      Ym-1[h-h   ] - Ym-1[h-h ]
                                   k+m            k
                   Y (h) = -------------------------
                    m           h  - h
                                 k    k+m

        and when h -> 0:

                         k         k+m     k+1       k
                    k   Y   [0-h0/2   ] - Y   [0-h0/2 ]
                         m-1               m-1
                   Y  = -------------------------------
                    m             k       k+m
                              h0/2  - h0/2

        Dividing by h0/2, we obtain:

                                 m  k+1    k
                           k    2  Ym-1 - Ym-1
                          Y  =  -------------
                           m       m
                                  2  - 1
                                            k
          Observing how the table elements Ym are calculated, we can see
        that it is not necessary to keep in memory a table with two dimen-
                                   k-1      k-2           0             k
        sions, we first calculate Y1, then Y2  ... until Yk. The table Ym
        is then replaced by the one-dimension table, Tk:

                                        k
           iteration #     step        Ym             Tk
           ---------------------------------------------------
                                      0          
               0           h         Y0            T1
               
                                      1  0
               1           h/2       Y0 Y1         T2  T1

                                      2  1  0
               2           h/4       Y0 Y1 Y2      T3  T2 T1

             ----         ----       ---------     ---------
                              k       k      0
               k           h/2       Y0 ... Yk     Tk ... T1

           ---------------------------------------------------

          To fully understand the mechanism, we recommend the reader to
                                            0   1      0   2
        manually calculate a few elements: Y1, Y1 and Y2, Y1, etc.

        Or else:
                               m
                              2  T    - T
                                  K=1    k
                        T  = -------------    (k = L...1)
                         k       m
                                2  - 1

          To sum up, we have in the program EULROMB the following steps:

          1. Define the function f(x,y), the starting values, X(1), Y(1),
             the starting step H, the maximum error ER, the max. number
             of subdivisions LA and the number of calculations NC.

                             NC = (X    - X ) / H
                                    n+1    1

          2. Approximate Yn+1 by Euler's method:

                   0                               0
                  Y0 = Yn + h f(xn,yn)  and  T1 = Y0

          3. Calculate the Tk elements of the extrapolation table; we
             initialize the iteration with L=1.

             a) XC = X  ; YC = Y        
                      n         n

                approximate Y    by Euler's method successively dividing
                             n+1
                the step by 2   ;   T    is the value of Y    given by
                                     L+1                  n+1
                Euler.

             b) Initialize the column count m=1; k=L is the T index.

                           k
             c) calculate Ym or rather T  using:
                                        k
                         m            
                        2  T    - T
                            k+1    k
                   T = -------------
                    k      m
                          2  - 1

             d) if k<1 goto subroutine f.

             e) test convergence: if |Tk - Tk+1| < ER, the iteration is
                finished, otherwise increment m, decrement k, go to step c)

             f) if L<Lmax, increment L then goto step a) else display a
                message "No convergence!"

             g) if n < N, increment n, go to 2. else end program.


             Note: the program EULROMB is made such as the solution Yk is
                   stored for all the points (this is not always necessary).


             [From BIBLI 04].

SOURCE WEBSITE: http://jean-pierre.moreau.pagesperso-orange.fr/eqdiff.html

Tuesday, August 29, 2017

VB6 Support for PNG, TIF, GIF Animation

LaVolpe is literally a genius and a super-advanced programmer who has revolutionized VB6 graphics. It revolutionized the VB6 graphics so much that Microsoft with all their army pay close attention to his methods and style. Yes, that means VB6 always had (and has) better graphics than Microsoft's newer products. LaVolpe is one of the people to whom we (The VB6 community) owe our superiority over Microsoft in recent years.

Download from me

Download from source



Status: Updated 14 Dec 2015. Fix rendering bugs discussed in post #37
Updated 11 Nov 2015 to include support for non-placeable WMFs

We all know VB is old and doesn't support many of the common image formats, specifically: TIF, PNG and alpha-blended icons. The attached class enables VB to support those formats and a bit more. There is no usercontrol involved. Just a bit of subclassing (IDE safe) to enable the images to be rendered with APIs that do support the alpha channel: AlphaBlend and DrawIconEx. The class is not limited to just Image controls, can be applied to most (if not all) of VB's picture, icon, and other image properties.

Image formats supported. The 'includes' below are in addition to what VB supports:
:: BMP. Includes 32bpp alpha and/or premultiplied. Includes those created with v4 & v5 of the bitmap info header
:: GIF. Includes methods to navigate frames
:: JPG. Includes CMYK color space
:: ICO,CUR. Includes 32bpp alpha and/or PNG-encoded Vista-type icons
:: WMF. Includes non-placeable
:: EMF
:: PNG
:: TIF. Both single page & multi-page. Supported compression schemes depend on version of GDI+ installed
See post #8 below for tips regarding design-view usage

Important to remember. This is still a VB stdPicture object. All of its properties and methods are supported, including Width, Height, Type, Render, Handle, etc.





The class methods/properties can be categorized as follows:

I. VB Replacement Methods
- LoadPicture: Unicode supported. Has options to keep original image format when loading. This enables you to navigate animated GIF frames and multipage TIFs. Cached data allows you to save that original data to file
- SavePicture: Unicode supported. Has options to always save to bitmap or save original format if it exists
- PictureType: Returns the original image format, i.e., GIF, PNG, TIF, etc, if that format exists
note: LoadPicture & SavePicture both accept byte array as source/destination medium

II. Management Methods
- IsManaged: Return whether the class is rendering the image or VB is
- UnManage: Simply unsubclasses a managed image
- HasOriginalFormat: Return whether or not any Picture is caching original image format data

III. Navigation Methods
- SubImageCount: Returns the number of frames/pages contained within the managed image
- SubImageIndex: Returns the current frame/page displayed by the managed image
- SubImage: Returns a stdPicture object of the requested frame/page
- GetGIFAnimationInfo: Returns an array containing loop count and each frame's display interval for animated GIFs

Quickly. How does this class do it?

1. It creates thunks to subclass the VB picture objects. The source code of the thunks are provided in the package. These thunks draw the image using the APIs AlphaBlend or DrawIconEx. For you ASM gurus out there: I'm a noob with ASM. I'm sure others can rewrite that ASM code to be more efficient.

2. To support AlphaBlend API, images may be converted to a 32bit bitmap having the pixel data premultiplied against the alpha channel. These bitmaps will never display correctly outside of the class. Therefore the class' SavePicture function allows you to create a copy of the bitmap that can be displayed outside of the class. This copy can be passed to the clipboard and/or drag/drop methods of your project.

3. GDI+ is relied upon to parse/load TIF and PNG. It is also used to support JPG in CMYK format, non-placeable WMFs and multi-frame GIF rendering. GDI+ is a requirement for the class. If it fails to load or thunks fail to be created, the class will silently fall back to standard VB methods and VB restrictions.

The transparency displayed by the image control is not faked. It is true transparency and the attached test project should easily prove that. Important to understand. TIF and PNG support is not available at design-time. This is because the class code isn't activated until run-time. Some motivated individuals out there could easily create a windowless usercontrol that hosts an image control (and this class) that could support all formats at design-time. Just a thought and subtle prod.

The class can be expanded by those willing to put in the effort. Ideas would be to incorporate GDI+ high quality scaling, conversion from one image format to another, image effects like rotation, blur, and more. Other image formats can easily be supported from outside the class. If you can parse/render that new format to a 32bpp bitmap, then you can use the class' LoadPicture to display that image. Have fun.
Tip: If modifying the thunks or code, recommend identifying these new versions using a different key. Throughout the code, you can find a couple instances of this text: IPIC+Thunker. Change those instances to reflect something else, i.e., "IPIC+Thunker.v2" so you can distinguish between your versions.

When compiled, VB can behave differently vs when uncompiled. Some differences are subtle, others are not. Here's one that is key for animating GIFs. In the test project posted below, the animation works because VB caches the GIF format for the GIF that was loaded into Image1 during design-time. During run-time that info is still cached by VB so the class can extract the entire GIF. But when you compile the project, the GIF no longer animates. Why? Well, when compiled, the GIF information is lost. VB no longer caches it. This can be proven in a simple test project. Add an image control and button. Add a GIF or JPG to that image control. Add the following code behind the button. Click the button when project is running compiled and uncompiled. Different results. The workaround is simply to save GIFs you want to animate in a resource file and load the GIF from that.

Known Limitations:
1) VB-like scaling is in play. Did not spend the time to incorporate high-quality GDI+ scaling
2) Do not send any Picture object to the class if it was retrieved via Clipboard.GetData. Why? Well, on 32 bit screen displays, the alpha channel can be filled with garbage. This garbage will likely be misinterpreted as transparency information
3) If expecting to use the .Render method of the picture object (icons only), do not keep the icon picture object returned from this class as Managed. After loading the icon, call UnManage. See Post #37 for more

Source:VBForums