meschach tutorial.htm

来自「C语言版本的矩阵库」· HTM 代码 · 共 870 行 · 第 1/3 页

    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);
      ......
}
</PRE>Normally, these registered workspace variables remain allocated. However, 
to implement the ``deallocate on exit'' approach, use the following code: <PRE>  ......
mem_stat_mark(1);
rk4(...,x,...)
mem_stat_free(1);
  ......
</PRE>To keep the workspace vectors allocated for the duration of a loop, but 
then deallocated, use <PRE>  ......
mem_stat_mark(1);
for (i = 0; i &lt; N; i++ )
    rk4(...,x,...);
mem_stat_free(1);
  ......
</PRE>The number used in the <TT>mem_stat_mark()</TT> and 
<TT>mem_stat_free()</TT> calls is the <I>workspace group number</I>. The call 
<TT>mem_stat_mark(1);</TT> designates 1 as the current workspace group number; 
the call <TT>mem_stat_free(1);</TT> deallocates (and sets to NULL) all static 
workspace variables registered as belonging to workspace group 1. 
<H2><A name=vec_ops>Simple vector operations: An RK4 routine</A></H2>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
<P><I>x'=f(t,x), x(t_0)=x_0</I>
<P>for <I>x(t_i)</I>, <I>i=1,2,3,...</I> where <I>t_i=t_0+i h</I> and <I>h</I> 
is the step size. The formulae for the 4th order Runge--Kutta method are: 
<UL>
  <LI>x_{i+1}=x_i+ (h/6){v_1+2v_2+2v_3+v_4} where 
  <LI>v_1=f(t_i,x_i) 
  <LI>v_2=f(t_i+ h/2,x_i+(h/2) v_1) 
  <LI>v_3=f(t_i+ h/2,x_i+(h/2) v_2) 
  <LI>v_4=f(t_i+h,x_i+h v_3) </LI></UL>where the <I>v_i</I> are vectors. The 
procedure for implementing this method (<TT>rk4()</TT>) will be passed (a 
pointer to) the function <I>f</I>; the implementation of <I>f</I> 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 <I>f</I> will also be passed an already 
allocated vector to be filled in with the appropriate values. The procedure 
<TT>rk4()</TT> will also be passed the current time <I>t</I>, the step size 
<I>h</I>, and the current value for <I>x</I>. The time after the step will be 
returned by <TT>rk4()</TT>. The code that does this follows. <PRE>#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, v2, etc. are of the correct size */
    v1   = v_resize(v1,x-&gt;dim);
    v2   = v_resize(v2,x-&gt;dim);
    v3   = v_resize(v3,x-&gt;dim);
    v4   = v_resize(v4,x-&gt;dim);
    temp = v_resize(temp,x-&gt;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 */
}
</PRE>Note that the last parameter of <TT>f()</TT> is where the <I>output</I> is 
placed. Often this can be <TT>NULL</TT> 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 <I>v_1,...,v_4</I> are created to have 
the correct size in the lines <PRE>v1 = v_resize(v1,x-&gt;dim);
v2 = v_resize(v2,x-&gt;dim);
....
</PRE>Here <TT>v_resize(v,dim)</TT> resizes the <TT>VEC</TT> structure 
<TT>v</TT> to hold a vector of length <TT>dim</TT>. If <TT>v</TT> is initially 
NULL, then this creates a new vector of dimension <TT>dim</TT>, just as 
<TT>v_get(dim)</TT> would do. For the above piece of code to work correctly, 
<TT>v1</TT>, <TT>v2</TT> etc., must be initialised to be NULL vectors. This is 
done by the declaration <PRE>static VEC *v1=VNULL, *v2=VNULL, *v3=VNULL, *v4=VNULL;
</PRE>or <PRE>static VEC *v1, *v2, *v3, *v4;
</PRE>The operations of vector addition and scalar addition are really the only 
<I>vector</I> operations that need to be performed in <TT>rk4</TT>. Vector 
addition is done by <PRE>v_add(v1,v2, out)</PRE>, where <TT>out=v1+v2</TT>, and scalar 
multiplication by <TT>sv_mlt(scale,v,out)</TT>, where <TT>out=scale*v</TT>. 
These can be combined into a single operation 
<TT>v_mltadd(v1,v2,scale,out)</TT>, where <TT>out=v1+scale*v2</TT>. As many 
operations in numerical mathematics involve accumulating scalar multiples, this 
is an extremely useful operation, as we can see above. For example: <PRE>v_mltadd(x,v1,0.5*h,temp);    /* temp = x+.5*h*v1 */
</PRE>We also need a number of ``utility'' operations. For example 
<TT>v_copy(in, out)</TT> copies the vector <TT>in</TT> to <TT>out</TT>. There is 
also <TT>v_zero(v)</TT> to zero a vector <TT>v</TT>. 
<P>Here is an implementation of the function <I>f</I> for simple harmonic 
motion: <PRE>/* 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-&gt;dim != 2 || out-&gt;dim != 2 )
        error(E_SIZES,"f");

    out-&gt;ve[0] = x-&gt;ve[1];
    out-&gt;ve[1] = - x-&gt;ve[0];

    return out;
}
</PRE>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 
<TT>f()</TT> 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 <TT>i</TT>th component of a vector <TT>x</TT> is 
<TT>x-&gt;ve[i]</TT>, and indexing is zero-relative (i.e., the ``first'' 
component is component 0). The ODE described above is for simple harmonic 
motion: <I>x_0'=x_1</I>, <I>x_1'=-x_0</I>, or equivalently, <I>x_0''+x_0=0</I>. 
Here is the main program: <PRE>#include <STDIO.H>
#include "matrix.h"

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

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

    printf("# At time 
    v_output(x);
    while ( t &lt; t_fin )
    {
        t = rk4(f,t,x,min(h,t_fin-t));/* new t is returned */
        printf("# At time 
        v_output(x);
        t += h;
    }
}
</PRE>Here the initial values are entered as a vector by <TT>v_input()</TT>. If 
<TT>v_input()</TT> is passed a vector, then this vector will be used to store 
the input, and this vector has the size that <TT>x</TT> had on entry to 
<TT>v_input()</TT>. The original values of <TT>x</TT> are also used as a prompt 
on input from a tty. If a <TT>NULL</TT> is passed to <TT>v_input()</TT> then 
<TT>v_input()</TT> 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 <TT>f()</TT> is expecting) we use <PRE>x = v_get(2);     x = v_input(x);
</PRE>
<P>To compile the program under Unix<I>(TM)</I>, if it is in a file 
<TT>tutorial.c</TT> is: <PRE>cc -o tutorial tutorial.c meschach.a
</PRE>or, if you have an ANSI compiler, <PRE>cc -DANSI_C -o tutorial tutorial.c meschach.a
</PRE>Here is a sample session with the above program: <PRE>
  ......
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 
</PRE>By way of comparison, the state at <I>t=1</I> for the true solution is 
<I>x_0(1)=0.5403023058</I>, <I>x_1(1)=-0.8414709848</I>. The ``<TT>b</TT>'' that 
is typed in entering the <TT>x</TT> 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 
``<TT>f</TT>'' for skipping over the vector's components, which are then 
unchanged. If an incorrectly sized initial value vector <TT>x</TT> is given, the 
error handler comes into action: <PRE>
  ......
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()
Sorry, aborting program

</PRE>The error handler prints out the error message giving the source code file 
and line number as well as the function name where the error was raised. The 
relevant section of <TT>f()</TT> in file <TT>test1.c</TT> is: <PRE>if ( x-&gt;dim != 2 || out-&gt;dim != 2 )
    error(E_SIZES,"f");               /* line 79 */
</PRE>The standard routines in this system perform error checking of this type, 
and also checking for undefined results such as division by zero in the routines 
for solving systems of linear equations. There are also error messages for 
incorrectly formatted input and end-of-file conditions. 
<P>To round off the discussion of this program, note that we have seen 
interactive input of vectors. If the input file or stream is not a tty (e.g., a 
file, a pipeline or a device) then it expects the input to <I>have the same form 
as the output for each of the data structures</I>. Each of the input routines 
(<TT>v_input()</TT>, <TT>m_input()</TT>, <TT>px_input()</TT>) skips over 
``comments'' in the input data, as do the macros <TT>input()</TT> and 
<TT>finput()</TT>. Anything from a `#' to the end of the line (or EOF) is 
considered to be a comment. For example, the initial value problem could be set 
up in a file <TT>ivp.dat</TT> as: <PRE># Initial time
0
# Final time
1
# Solution is x(t) = (cos(t),-sin(t))
# x(0) =
Vector: dim: 2
1       0
# Step size
0.1
</PRE>The output of the above program with the above input (from a file) gives 
essentially the same output as shown above, except that no prompts are sent to 
the screen. 
<H2><A name=list_args>Using routines for lists of arguments</A></H2>Some of the 
most common routines have vaariants that take a variable number of arguments. 
These are the routines <TT>..get_vars()</TT>, <TT>.._resize_vars()</TT> and 
<TT>.._free_vars()</TT>. These correspond to the the basic routines 
<TT>.._get()</TT>, <TT>.._resize()</TT> and <TT>.._free()</TT> respectively. 
Also there is the <TT>mem_stat_reg_vars()</TT> routine which registers a list of 
static workspace variables; this corresponds to <TT>mem_stat_reg_list()</TT> for 
a single variable. Here is an example of how to use these functions. This 
example, also uses the routine <TT>v_linlist()</TT> to compute a linear 
combinartion. Note that the code is much more compact, but don't forget that 
these ``<TT>..._vars()</TT>'' routines usually need the address-of operator 
``<TT>&amp;</TT>'' and NULL termination of the arguments for these to work 
correctly. <PRE>#include "matrix.h"

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

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

    /* ensure that v1, v2 etc. are of the correct size */
    v_resize_vars(x-&gt;dim,&amp;v1,&amp;v2,&amp;v3,&amp;v4,&amp;temp,NULL);
    /* register workspace variables */
    mem_stat_reg_vars(0,TYPE_VEC,&amp;v1,&amp;v2,&amp;v3,&amp;v4,&amp;temp,NULL);
    /* end of memory allocation */
    (*f)(t,x,v1);             v_mltadd(x,v1,0.5*h,temp);
    (*f)(t+0.5*h,temp,v2);    v_mltadd(x,v2,0.5*h,temp);
    (*f)(t+0.5*h,temp,v3);    v_mltadd(x,v3,h,temp);
    (*f)(t+h,temp,v4);

    /* now add: temp = v1+2*v2+2*v3+v4 */