tutorial.txt

C语言版本的矩阵库

  rk4(...,x,...)
  {
     static VEC *v1=VNULL, *v2=VNULL, *v3=VNULL, *v4=VNULL, *temp=VNULL;
     .......
     v1   = v_resize(v1,x->dim);
     v2   = v_resize(v2,x->dim);
     v3   = v_resize(v3,x->dim);
     v4   = v_resize(v4,x->dim);
     temp = v_resize(temp,x->dim);
     .......
  }

   The intention is that the rk4() routine is called repeatedly with the
same size x vector.  It then doesn't make as much sense to allocate v1, v2
etc.  whenever the function is called.  Instead, v_resize() only performs
memory allocation if the memory already allocated to v1, v2 etc. is smaller
than x->dim.

   The vectors v1, v2 etc. are declared to be static to ensure that their
values are not lost between function calls.  Variables that are declared
static are set to NULL or zero by default.  So the declaration of v1, v2,
etc., could be

  static VEC *v1, *v2, *v3, *v4, *temp;

   This strategy of resizing static workspace variables is not so useful if
the object being allocated is extremely large.  The previous ``allocate and
deallocate'' strategy is much more efficient for memory in those
circumstances.  However, the following section shows how to get the best of
both worlds.


2.4  REGISTRATION OF WORKSPACE

   From version 1.2 onwards, workspace variables can be registered so that
the memory they reference can be freed up on demand.  To do this, the
function containing the static workspace variables has to include calls to
MEM_STAT_REG(var,type) where var is a pointer to a Meschach data type (such
as VEC or MAT).  This call should be placed after the call to the
appropriate resize function.  The type parameter should be a TYPE_... macro
where the ``...'' is the name of a Meschach type such as VEC or MAT.  For
example,

  rk4(...,x,...)
  {
     static VEC *v1, *v2, *v3, *v4, *temp;
       .......
     v1   = v_resize(v1,x->dim);
     MEM_STAT_REG(v1,TYPE_VEC);
     v2   = v_resize(v2,x->dim);
     MEM_STAT_REG(v2,TYPE_VEC);
       ......
  }

Normally, these registered workspace variables remain allocated.  However,
to implement the ``deallocate on exit'' approach, use the following code:

  ......
  mem_stat_mark(1);
  rk4(...,x,...)
  mem_stat_free(1);
  ......

   To keep the workspace vectors allocated for the duration of a loop, but
then deallocated, use

  ......
  mem_stat_mark(1);
  for (i = 0; i < N; i++ )
    rk4(...,x,...);
  mem_stat_free(1);
  ......

The number used in the mem_stat_mark() and mem_stat_free() calls is the
workspace group number.  The call mem_stat_mark(1) designates 1 as the
current workspace group number; the call mem_stat_free(1) deallocates (and
sets to NULL) all static workspace variables registered as belonging to
workspace group 1.



3.  SIMPLE VECTOR OPERATIONS: AN RK4 ROUTINE

   The main purpose of this example is to show how to deal with vectors and
to compute linear combinations.

   The problem here is to implement the standard 4th order Runge-Kutta
method for the ODE

  x'=f(t,x), x(t_0)=x_0 

for x(t_i), i=1,2,3, where t_i=t_0+i*h and h is the step size.

   The formulae for the 4th order Runge-Kutta method are:

	x_i+1 = x_i+ h/6*(v_1+2*v_2+2*v_3+v_4),
where
	v_1 = f(t_i,x_i)
	v_2 = f(t_i+h, x_i+h*v_1)
	v_3 = f(t_i+h, x_i+h*v_2)
	v_4 = f(t_i+h, x_i+h*v_3)

where the v_i are vectors.

   The procedure for implementing this method (rk4()) will be passed (a
pointer to) the function f. The implementation of f could, in this system,
create a vector to hold the return value each time it is called.  However,
such a scheme is memory intensive and the calls to the memory allocation
functions could easily dominate the time performed doing numerical
computations.  So, the implementation of f will also be passed an already
allocated vector to be filled in with the appropriate values.

   The procedure rk4() will also be passed the current time t, the step
size h, and the current value for x.  The time after the step will be
returned by rk4().

The code that does this follows.


  #include "matrix.h"

  /* rk4 - 4th order Runge-Kutta method */
  double rk4(f,t,x,h)
  double t, h;
  VEC    *(*f)(), *x;
  {
     static VEC *v1=VNULL, *v2=VNULL, *v3=VNULL, *v4=VNULL;
     static VEC *temp=VNULL;

     /* do not work with NULL initial vector */
     if ( x == VNULL )
        error(E_NULL,"rk4");

     /* ensure that v1, ..., v4, temp are of the correct size */
     v1   = v_resize(v1,x->dim);
     v2   = v_resize(v2,x->dim);
     v3   = v_resize(v3,x->dim);
     v4   = v_resize(v4,x->dim);
     temp = v_resize(temp,x->dim);

     /* register workspace variables */
     MEM_STAT_REG(v1,TYPE_VEC);
     MEM_STAT_REG(v2,TYPE_VEC);
     MEM_STAT_REG(v3,TYPE_VEC);
     MEM_STAT_REG(v4,TYPE_VEC);
     MEM_STAT_REG(temp,TYPE_VEC);
     /* end of memory allocation */

     (*f)(t,x,v1);         /* most compilers allow: f(t,x,v1); */
     v_mltadd(x,v1,0.5*h,temp);    /* temp = x+.5*h*v1 */
     (*f)(t+0.5*h,temp,v2);
     v_mltadd(x,v2,0.5*h,temp);    /* temp = x+.5*h*v2 */
     (*f)(t+0.5*h,temp,v3);
     v_mltadd(x,v3,h,temp);        /* temp = x+h*v3 */
     (*f)(t+h,temp,v4);

     /* now add: v1+2*v2+2*v3+v4 */
     v_copy(v1,temp);              /* temp = v1 */
     v_mltadd(temp,v2,2.0,temp);   /* temp = v1+2*v2 */
     v_mltadd(temp,v3,2.0,temp);   /* temp = v1+2*v2+2*v3 */
     v_add(temp,v4,temp);          /* temp = v1+2*v2+2*v3+v4 */

     /* adjust x */
     v_mltadd(x,temp,h/6.0,x);     /* x = x+(h/6)*temp */

     return t+h;                   /* return the new time */
  }


   Note that the last parameter of f() is where the output is placed.
Often this can be NULL in which case the appropriate data structure is
allocated and initialised.  Note also that this routine can be used for
problems of arbitrary size, and the dimension of the problem is determined
directly from the data given.  The vectors v_1,...,v_4 are created to have
the correct size in the lines

  ....
  v1 = v_resize(v1,x->dim);
  v2 = v_resize(v2,x->dim);
  ....

   Here v_resize(v,dim) resizes the VEC structure v to hold a vector of
length dim.  If v is initially NULL, then this creates a new vector of
dimension dim, just as v_get(dim) would do.  For the above piece of code to
work correctly, v1, v2 etc., must be initialised to be NULL vectors.  This
is done by the declaration

  static VEC *v1=VNULL, *v2=VNULL, *v3=VNULL, *v4=VNULL;

or

  static VEC *v1, *v2, *v3, *v4;

The operations of vector addition and scalar addition are really the only
vector operations that need to be performed in rk4.  Vector addition is
done by v_add(v1,v2,out), where out=v1+v2, and scalar multiplication by
sv_mlt(scale,v,out), where out=scale*v.

These can be combined into a single operation v_mltadd(v1,v2,scale,out),
where out=v1+scale*v2.  As many operations in numerical mathematics involve
accumulating scalar multiples, this is an extremely useful operation, as we
can see above.  For example:

  v_mltadd(x,v1,0.5*h,temp);    /* temp = x+0.5*h*v1 */

   We also need a number of ``utility'' operations.  For example v_copy(in,
out) copies the vector in to out.  There is also v_zero(v) to zero a vector
v.

   Here is an implementation of the function f for simple harmonic motion:

  /* f - right-hand side of ODE solver */
  VEC	*f(t,x,out)
  VEC	*x, *out;
  double	t;
  {
    if ( x == VNULL || out == VNULL )
        error(E_NULL,"f");
    if ( x->dim != 2 || out->dim != 2 )
        error(E_SIZES,"f");

    out->ve[0] = x->ve[1];
    out->ve[1] = - x->ve[0];

    return out;
  }

  As can be seen, most of this code is error checking code, which, of
course, makes the routine safer but a little slower.  For a procedure like
f() it is probably not necessary, although then the main program would have
to perform checking to ensure that the vectors involved have the correct
size etc.  The ith component of a vector x is x->ve[i], and indexing is
zero-relative (i.e., the ``first'' component is component 0).  The ODE
described above is for simple harmonic motion:
 
	x_0'=x_1 ,  x_1'=-x_0 , or equivalently,  x_0''+ x_0 = 0 .

  Here is the main program:


  #include <stdio.h>
  #include "matrix.h"

  main()
  {
    VEC        *x;
    VEC        *f();
    double     h, t, t_fin;
    double     rk4();

    input("Input initial time: ", "%lf", &t);
    input("Input final time: ",  "%lf", &t_fin);
    x = v_get(2);        /* this is the size needed by f() */
    prompter("Input initial state:\n");	x = v_input(VNULL);
    input("Input step size: ", "%lf", &h);

    printf("# At time %g, the state is\n",t); 
    v_output(x);
    while ( t < t_fin )
    {
        t = rk4(f,t,x,min(h,t_fin-t));   /* new t is returned */
        printf("# At time %g, the state is\n",t);
        v_output(x);
	t += h;
    }
  }

   The initial values are entered as a vector by v_input().  If v_input()
is passed a vector, then this vector will be used to store the input, and
this vector has the size that x had on entry to v_input().  The original
values of x are also used as a prompt on input from a tty.  If a NULL is
passed to v_input() then v_input() will return a vector of whatever size
the user inputs.  So, to ensure that only a two-dimensional vector is used
for the initial conditions (which is what f() is expecting) we use

	x = v_get(2);     x = v_input(x);

   To compile the program under Unix, if it is in a file tutorial.c:

	cc -o tutorial tutorial.c meschach.a

or, if you have an ANSI compiler,

	cc -DANSI_C -o tutorial tutorial.c meschach.a

   Here is a sample session with the above program: 

 tutorial

  Input initial time: 0
  Input final time: 1
  Input initial state:
  Vector: dim: 2
  entry 0: -1
  entry 1: b
  entry 0: old             -1 new: 1
  entry 1: old              0 new: 0
  Input step size: 0.1
  At time 0, the state is
  Vector: dim: 2
             1              0 
  At time 0.1, the state is
  Vector: dim: 2
    0.995004167  -0.0998333333 
      .................
  At time 1, the state is
  Vector: dim: 2
    0.540302967   -0.841470478 

   By way of comparison, the state at t=1 for the true solution is
	 x_0(1)=0.5403023058 , x_1(1)=-0.8414709848 .  
The ``b'' that is typed in entering the x vector allows the user to alter
previously entered components. In this case once this is done, the user is
prompted with the old values when entering the new values.  The user can
also type in ``f'' for skipping over the vector's components, which are
then unchanged.  If an incorrectly sized initial value vector x is given,
the error handler comes into action:

  Input initial time: 0
  Input final time: 1
  Input initial state:
  Vector: dim: 3
  entry 0: 3
  entry 1: 2
  entry 2: -1
  Input step size: 0.1
  At time 0, the state is
  Vector: dim: 3
             3              2             -1 

  "tutorial.c", line 79: sizes of objects don't match in function f()