Anubis / Anubis

   
 *Project*                             The Anubis Project
   
 *Title*               Drawing solutions of first order differential equations
                            using the fourth order Runge-Kutta method. 
   
 *Copyright*                     Copyright (c) Alain Prout� 2005. 

 *Released*
   
 *Author*       Alain Prout�
   

 *Overview* 
   The program defined  in this file opens a  host window.  When a point of  the window is
   clicked upon,  the solution of a  given first order differential  equation which passes
   through this point is drawn.
   
   First order differential equations (when solved  with respect to y') have the following
   form:
   
                                            y' = f(x,y)
   
   where y is the  unknown function, x the free (time) variable,  and f a function defined
   on a  subset of  R^2. A  solution of  this equation is  a function  g (defined  on some
   interval of the reals), such that for all x in this interval, we have:
   
                                         g'(x) = f(x,g(x))
   
   The Cauchy-Lipschitz theorem asserts that if  f is locally lipschitzian relative to its
   second variable, then,  by every point of the  domain of f, passes one  and exactly one
   maximal (i.e. non extendable) solution of the equation. 
   
   Below is the function f used in this example (you may customize). 

read tools/basis.anubis   
read tools/maybefloat.anubis

   
   
   
define Maybe(Float)
   f
     (
       Maybe(Float)  x, 
       Maybe(Float)  y
     ) =   
   x - y*y.
   sin(x*y). 
   _3 * ((y*y)^(_1/_3)).
   tan(y). 
   x*x - y*y.  
   x - y*x*y.
   -y*y*y. 
   abs(x)^y. 
   sin(y). 
   _1/y. 

   

   
   x+y*y*y. 
   x*x + y. 
   x*x - y*y. 

   
   
   
  
   The Runge-Kutta  method enables the  computation of points  on the graph of  a solution
   starting at some point  of this graph.  Starting at a point  (x_0,y_0) in the domain of
   f,  we compute the  points (x_n,y_n)  belonging to  the graph  of the  solution passing
   through (x_0,y_0) as follows:
   
         x_(n+1) = x_n + h
         k_1 = hf(x_n,y_n)
         k_2 = hf(x_n + h/2, y_n + k_1/2)
         k_3 = hf(x_n + h/2, y_n + k_2/2)
         k_4 = hf(x_n + h, y_n + k_3)
         y_(n+1) = y_n + k_1/6 + k_2/3 + k_3/3 + k_4/6
   
   The error  made by this method  against the actual solution  is of the order  of h^5 at
   each step.
   
 define Float
   to_Float
     (
       Word32 x
     ) =
   to_Float(to_Int(x)). 
   
   
define Float
   _to_Float
     (
       Int x
     ) =
   to_Float(truncate_to_Word32(x)). 
   
   The next function computes y_(n+1) from x_n, y_n and h, using f. 
   
define Maybe(Float)
   runge_kutta_4
     (
       Maybe(Float)                                          x_n,
       Maybe(Float)                                          y_n,
       Maybe(Float)                                          h, 
       (Maybe(Float) x, Maybe(Float) y) -> Maybe(Float)      f
     ) =
   with h2 = h/success(2), 
        k_1 = h * f(x_n,y_n),
        k_2 = h * f(x_n + h2, y_n + k_1/success(2)), 
        k_3 = h * f(x_n + h2, y_n + k_2/success(2)),
        k_4 = h * f(x_n + h, y_n + k_3),
     y_n + ((k_1+k_4)/success(6)) + ((k_2+k_3)/success(3)). 

   
   A point of  coordinates (x,y) must be  painted into the buffer (of  type HostImage). We
   have two scaling numbers 'h_scale' and 'v_scale'. 
   
define RGB white = rgb(255,255,255).    
   

   
define Bool
   paint
     (
       HostImage      buffer,
       Maybe(Float)   x, 
       Maybe(Float)   y, 
       Maybe(Float)   h_scale, 
       Maybe(Float)   v_scale,
       Int            xmax, 
       Int            ymax
     ) =
   with x1 = x/h_scale, y1 = y/v_scale, 
   if x1 is success(x2) 
   then if y1 is success(y2)
        then with x3 = integral_part(x2), 
                  y3 = integral_part(y2), 
               if (0 =< x3 & x3 < xmax & 0 =< y3 & y3 < ymax)
               then (paint_rectangle(buffer,rect(x3,y3,x3+1,y3+1),white); true)
               else false
        else false
   else false.
   

   Now, if we want to draw a solution we must iterate the above function.
   
define One
   draw_solution 
     (
       Maybe(Float)                                          x, 
       Maybe(Float)                                          y,
       Maybe(Float)                                          h, 
       (Maybe(Float) x, Maybe(Float) y) -> Maybe(Float)      f,
       Maybe(Float)                                          h_scale,
       Maybe(Float)                                          v_scale,
       Maybe(Float)                                          x_trans, 
       Maybe(Float)                                          y_trans, 
       HostImage                                             buffer,
       Int                                                   xmax,
       Int                                                   ymax
     ) =
   if y is 
     {
       failure then unique,    // overflow
       success(_) then 
         if paint(buffer,
                  x,
                  success(_to_Float(ymax))*v_scale - y,
                  h_scale,
                  v_scale,
                  xmax,
                  ymax)
         then draw_solution(x+h,
                            runge_kutta_4(x,y,h,
                              (Maybe(Float) _x, Maybe(Float) _y) |-> 
                                f(x+x_trans,y+y_trans)),
                            h,
                            f,
                            h_scale,v_scale,
                            x_trans,y_trans, 
                            buffer,
                            xmax,ymax)
         else unique
     }.
   
   
   
   The next function makes  the paint_method. It just copies the content  of the buffer to
   the window, but only within the given rectangle. 
   
define (HostWindow hw, List(Rectangle) rs) -> One
   make_paint_method
     (
       HostImage buffer,
       Int x_trans, 
       Int y_trans,
       Int xmax, 
       Int ymax,
       Int hs, 
       Int vs
     ) =
   with red = rgb(255,0,0), 
   (HostWindow hw, List(Rectangle) rs) |->
     forget(map((Rectangle r) |-> paint_image(hw,r,0,0,buffer),rs));   // copy buffer to window within r
     paint_rectangle(hw,rect(
                             0,
                             ymax+y_trans-1,
                             xmax,
                             ymax+y_trans
                            ),red);
     paint_rectangle(hw,rect(
                             -x_trans+hs,
                             ymax+y_trans-4,
                             -x_trans+hs+1,
                             ymax+y_trans+4
                            ),red); 
     paint_rectangle(hw,rect(
                             -x_trans,
                             0,
                             -x_trans+1,
                             ymax
                            ),red);
     paint_rectangle(hw,rect(
                             -x_trans-4,
                             ymax+y_trans-vs-1,
                             -x_trans+4,
                             ymax+y_trans-vs
                            ),red). 
   
   
   The event handler.
   
   
define One
   handle_keys
     (
       KeyboardState ks, 
       KeyboardKey kk, 
       Var(Maybe(Float)) h_scale, 
       Var(Maybe(Float)) v_scale
     ) =
   if kk is 
     {
       unknown       then unique, 
       char(c)       then if c = '+' then (h_scale <- (*h_scale) * success(0.5);
                                           v_scale <- (*v_scale) * success(0.5)) else
                          if c = '-' then (h_scale <- (*h_scale) * success(2);
                                           v_scale <- (*v_scale) * success(2)) else
                          unique, 
       unicode(a,b)  then unique, 
       special(s)    then unique
     }. 
   
   
define (HostWindow hw, HostWindowEvent(One) e) -> List(Rectangle)
   make_event_handler
     (
       Rectangle default,
       HostImage buffer,
       Var(Maybe(Float)) h_scale, 
       Var(Maybe(Float)) v_scale, 
       Int h, 
       Int x_trans, 
       Int y_trans, 
       Int xmax, 
       Int ymax
     ) =
   (HostWindow hw, HostWindowEvent(One) e) |-> 
     if e is 
       {
         quit                    then [],
         expose                  then [default]
         pointer_entering        then [],
         pointer_leaving         then [],
         key_down(ks,kk)         then handle_keys(ks,
                                                  kk,
                                                  h_scale,
                                                  v_scale); [default],
         mouse_move(ks,x,y)      then [],
         mouse_click(ks,mc,x,y)  then if mc is left_down
                                      then (
                                           draw_solution(   // to the right (if h > 0)
                                             success(to_Float(x))**h_scale,
                                             success(_to_Float(ymax-to_Int(y)))**v_scale,
                                             success(_to_Float(h))**h_scale,
                                             f, 
                                             *h_scale,
                                             *v_scale,
                                             success(_to_Float(x_trans))**h_scale,
                                             success(_to_Float(y_trans))**v_scale,
                                             buffer,
                                             xmax,ymax); 
                                           draw_solution(   // to the left (if h > 0)
                                             success(to_Float(x))**h_scale,
                                             success(_to_Float(ymax-to_Int(y)))**v_scale,
                                             success(_to_Float(-h))**h_scale,
                                             f, 
                                             *h_scale,
                                             *v_scale,
                                             success(_to_Float(x_trans))**h_scale,
                                             success(_to_Float(y_trans))**v_scale,
                                             buffer,
                                             xmax,ymax); 
                                           [default]
                                           )
                                      else [],
         mouse_wheel(ks,d,x,y)   then [],
         tick                    then [],
         repaint(l)              then [],
         specific(_)             then []
       }. 
   
   
   Finally,  here is  our  program.  It  first creates  the  buffer, then  opens the  host
   window. If the  host window has been opened,  it must be shown by  'show', otherwise it
   remains invisible.
   
global define One
   runge_kutta
     (
       List(String) args
     ) =
   with iw = (Int)800, 
        ih = (Int)600, 
        default = rect(0,0,iw,ih), 
        hs = (Int)60, 
        vs = (Int)60, 
        step = (Int)1, 
        x_trans = (Int)-(iw\2), 
        y_trans = (Int)-(ih\2), 
   with buffer = to_host_image(create_rgba_image(iw,ih,rgba(10,10,10,255)),1),
   if open_host_window(rect(100,50,100+iw,50+ih),
                       "4-th order Runge Kutta method (please, click into the window)",
                       managed(not_resizable),
                       make_paint_method(buffer,x_trans,y_trans,iw,ih,hs,vs),
                       make_event_handler(default,buffer,
                                          var(success(1.0)/success(_to_Float(hs))), 
                                          var(success(1.0)/success(_to_Float(vs))), 
                                          step,
                                          x_trans,y_trans,
                                          iw,ih),
                       (List(HostWindowEvent(One)) l) |-> l) is 
     {
       failure then print("Cannot open window.\n"),
       success(hw) then show(hw)
     }.