home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 22 gnu / 22-gnu.zip / mesch12a.zip / tutorial.c < prev    next >
C/C++ Source or Header  |  1994-01-17  |  8KB  |  321 lines
  1.  
  2. /**************************************************************************
  3. **
  4. ** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  5. **
  6. **                 Meschach Library
  7. ** 
  8. ** This Meschach Library is provided "as is" without any express 
  9. ** or implied warranty of any kind with respect to this software. 
  10. ** In particular the authors shall not be liable for any direct, 
  11. ** indirect, special, incidental or consequential damages arising 
  12. ** in any way from use of the software.
  13. ** 
  14. ** Everyone is granted permission to copy, modify and redistribute this
  15. ** Meschach Library, provided:
  16. **  1.  All copies contain this copyright notice.
  17. **  2.  All modified copies shall carry a notice stating who
  18. **      made the last modification and the date of such modification.
  19. **  3.  No charge is made for this software or works derived from it.  
  20. **      This clause shall not be construed as constraining other software
  21. **      distributed on the same medium as this software, nor is a
  22. **      distribution fee considered a charge.
  23. **
  24. ***************************************************************************/
  25.  
  26. /* tutorial.c 10/12/1993 */
  27.  
  28. /* routines from Chapter 1 of Meschach */
  29.  
  30. static char rcsid[] = "$Id: tutorial.c,v 1.3 1994/01/16 22:53:09 des Exp $";
  31.  
  32. #include "matrix.h"
  33.  
  34. /* rk4 -- 4th order Runge--Kutta method */
  35. double rk4(f,t,x,h)
  36. double t, h;
  37. VEC    *(*f)(), *x;
  38. {
  39.    static VEC *v1=VNULL, *v2=VNULL, *v3=VNULL, *v4=VNULL;
  40.    static VEC *temp=VNULL;
  41.    
  42.    /* do not work with NULL initial vector */
  43.    if ( x == VNULL )
  44.      error(E_NULL,"rk4");
  45.  
  46.    /* ensure that v1, ..., v4, temp are of the correct size */
  47.    v1   = v_resize(v1,x->dim);
  48.    v2   = v_resize(v2,x->dim);
  49.    v3   = v_resize(v3,x->dim);
  50.    v4   = v_resize(v4,x->dim);
  51.    temp = v_resize(temp,x->dim);
  52.  
  53.    /* register workspace variables */
  54.    MEM_STAT_REG(v1,TYPE_VEC);
  55.    MEM_STAT_REG(v2,TYPE_VEC);
  56.    MEM_STAT_REG(v3,TYPE_VEC);
  57.    MEM_STAT_REG(v4,TYPE_VEC);
  58.    MEM_STAT_REG(temp,TYPE_VEC);
  59.    /* end of memory allocation */
  60.  
  61.    (*f)(t,x,v1); /* most compilers allow: "f(t,x,v1);" */
  62.    v_mltadd(x,v1,0.5*h,temp);    /* temp = x+.5*h*v1 */
  63.    (*f)(t+0.5*h,temp,v2);
  64.    v_mltadd(x,v2,0.5*h,temp);    /* temp = x+.5*h*v2 */
  65.    (*f)(t+0.5*h,temp,v3);
  66.    v_mltadd(x,v3,h,temp);        /* temp = x+h*v3 */
  67.    (*f)(t+h,temp,v4);
  68.    
  69.    /* now add: v1+2*v2+2*v3+v4 */
  70.    v_copy(v1,temp);              /* temp = v1 */
  71.    v_mltadd(temp,v2,2.0,temp);   /* temp = v1+2*v2 */
  72.    v_mltadd(temp,v3,2.0,temp);   /* temp = v1+2*v2+2*v3 */
  73.    v_add(temp,v4,temp);          /* temp = v1+2*v2+2*v3+v4 */
  74.    
  75.    /* adjust x */
  76.    v_mltadd(x,temp,h/6.0,x);     /* x = x+(h/6)*temp */
  77.    
  78.    return t+h;                   /* return the new time */
  79. }
  80.  
  81.  
  82.  
  83. /* rk4 -- 4th order Runge-Kutta method */
  84. /* another variant */
  85. double rk4_var(f,t,x,h)
  86. double t, h;
  87. VEC    *(*f)(), *x;
  88. {
  89.    static VEC *v1, *v2, *v3, *v4, *temp;
  90.    
  91.    /* do not work with NULL initial vector */
  92.    if ( x == VNULL )        error(E_NULL,"rk4");
  93.    
  94.    /* ensure that v1, ..., v4, temp are of the correct size */
  95.    v_resize_vars(x->dim, &v1, &v2, &v3, &v4, &temp, NULL);
  96.  
  97.    /* register workspace variables */
  98.    mem_stat_reg_vars(0, TYPE_VEC, &v1, &v2, &v3, &v4, &temp, NULL);
  99.    /* end of memory allocation */
  100.  
  101.    (*f)(t,x,v1);             v_mltadd(x,v1,0.5*h,temp);
  102.    (*f)(t+0.5*h,temp,v2);    v_mltadd(x,v2,0.5*h,temp);
  103.    (*f)(t+0.5*h,temp,v3);    v_mltadd(x,v3,h,temp);
  104.    (*f)(t+h,temp,v4);
  105.    
  106.    /* now add: temp = v1+2*v2+2*v3+v4 */
  107.    v_linlist(temp, v1, 1.0, v2, 2.0, v3, 2.0, v4, 1.0, VNULL);
  108.    /* adjust x */
  109.    v_mltadd(x,temp,h/6.0,x);     /* x = x+(h/6)*temp */
  110.    
  111.    return t+h;                   /* return the new time */
  112. }
  113.  
  114.  
  115. /* f -- right-hand side of ODE solver */
  116. VEC    *f(t,x,out)
  117. VEC    *x, *out;
  118. double    t;
  119. {
  120.    if ( x == VNULL || out == VNULL )
  121.      error(E_NULL,"f");
  122.    if ( x->dim != 2 || out->dim != 2 )
  123.      error(E_SIZES,"f");
  124.    
  125.    out->ve[0] = x->ve[1];
  126.    out->ve[1] = - x->ve[0];
  127.    
  128.    return out;
  129. }
  130.  
  131.  
  132. void tutor_rk4()
  133. {
  134.    VEC        *x;
  135.    VEC        *f();
  136.    double     h, t, t_fin;
  137.    double     rk4();
  138.    
  139.    input("Input initial time: ","%lf",&t);
  140.    input("Input final time: ",  "%lf",&t_fin);
  141.    x = v_get(2);    /* this is the size needed by f() */
  142.    prompter("Input initial state:\n");    x = v_input(VNULL);
  143.    input("Input step size: ",   "%lf",&h);
  144.    
  145.    printf("# At time %g, the state is\n",t);
  146.    v_output(x);
  147.    while (t < t_fin)
  148.    {
  149.       /* you can use t = rk4_var(f,t,x,min(h,t_fin-t)); */
  150.       t = rk4(f,t,x,min(h,t_fin-t));   /* new t is returned */
  151.       printf("# At time %g, the state is\n",t);
  152.       v_output(x);
  153.    }
  154. }
  155.  
  156.  
  157.  
  158.  
  159. #include "matrix2.h"
  160.  
  161. void tutor_ls()
  162. {
  163.    MAT *A, *QR;
  164.    VEC *b, *x, *diag;
  165.    
  166.    /* read in A matrix */
  167.    printf("Input A matrix:\n");
  168.    
  169.    A = m_input(MNULL);     /* A has whatever size is input */
  170.    
  171.    if ( A->m < A->n )
  172.    {
  173.       printf("Need m >= n to obtain least squares fit\n");
  174.       exit(0);
  175.    }
  176.    printf("# A =\n");       m_output(A);
  177.    diag = v_get(A->m);
  178.    /* QR is to be the QR factorisation of A */
  179.    QR = m_copy(A,MNULL);
  180.    QRfactor(QR,diag);   
  181.    /* read in b vector */
  182.    printf("Input b vector:\n");
  183.    b = v_get(A->m);
  184.    b = v_input(b);
  185.    printf("# b =\n");       v_output(b);
  186.    
  187.    /* solve for x */
  188.    x = QRsolve(QR,diag,b,VNULL);
  189.    printf("Vector of best fit parameters is\n");
  190.    v_output(x);
  191.    /* ... and work out norm of errors... */
  192.    printf("||A*x-b|| = %g\n",
  193.       v_norm2(v_sub(mv_mlt(A,x,VNULL),b,VNULL)));
  194. }
  195.  
  196.  
  197. #include <math.h>
  198. #include "iter.h"
  199.  
  200.  
  201. #define N 50
  202. #define VEC2MAT(v,m)  vm_move((v),0,(m),0,0,N,N);
  203.  
  204. #define PI 3.141592653589793116
  205. #define index(i,j) (N*((i)-1)+(j)-1)
  206.  
  207. /* right hand side function (for generating b) */
  208. double f1(x,y)
  209. double x,y;
  210. {
  211.   /* return 2.0*PI*PI*sin(PI*x)*sin(PI*y); */
  212.    return exp(x*y);
  213. }
  214.  
  215. /* discrete laplacian */
  216. SPMAT *laplacian(A)
  217. SPMAT *A;
  218. {
  219.    Real h;
  220.    int i,j;
  221.    
  222.    if (!A)
  223.      A = sp_get(N*N,N*N,5);
  224.  
  225.    for ( i = 1; i <= N; i++ )
  226.      for ( j = 1; j <= N; j++ )
  227.      {
  228.         if ( i < N )
  229.       sp_set_val(A,index(i,j),index(i+1,j),-1.0);
  230.         if ( i > 1 )
  231.       sp_set_val(A,index(i,j),index(i-1,j),-1.0);
  232.         if ( j < N )
  233.       sp_set_val(A,index(i,j),index(i,j+1),-1.0);
  234.         if ( j > 1 )
  235.       sp_set_val(A,index(i,j),index(i,j-1),-1.0);
  236.         sp_set_val(A,index(i,j),index(i,j),4.0);
  237.      }
  238.    return A;
  239. }
  240.  
  241. /* generating right hand side */
  242. VEC *rhs_lap(b)
  243. VEC *b;
  244. {
  245.    Real h,h2,x,y;
  246.    int i,j;
  247.    
  248.    if (!b)
  249.      b = v_get(N*N);
  250.  
  251.    h = 1.0/(N+1);      /* for a unit square */
  252.    h2 = h*h;
  253.    x = 0.0;
  254.    for ( i = 1; i <= N; i++ ) {
  255.       x += h;
  256.       y = 0.0;
  257.      for ( j = 1; j <= N; j++ ) {
  258.     y += h;
  259.     b->ve[index(i,j)] = h2*f1(x,y);
  260.      }
  261.    }
  262.    return b;
  263. }
  264.    
  265. void tut_lap()
  266. {
  267.    SPMAT *A, *LLT;
  268.    VEC *b, *out, *x;
  269.    MAT *B;
  270.    int num_steps;
  271.    FILE *fp;
  272.  
  273.    A = sp_get(N*N,N*N,5);
  274.    b = v_get(N*N);
  275.  
  276.    laplacian(A);
  277.    LLT = sp_copy(A);
  278.    spICHfactor(LLT);
  279.  
  280.    out = v_get(A->m);
  281.    x = v_get(A->m);
  282.  
  283.    rhs_lap(b);   /* new rhs */
  284.    iter_spcg(A,LLT,b,1e-6,out,1000,&num_steps);
  285.    printf("Number of iterations = %d\n",num_steps);
  286.  
  287.    /* save b as a MATLAB matrix */
  288.  
  289.    fp = fopen("laplace.mat","w");  /* b will be saved in laplace.mat */
  290.    if (fp == NULL) {
  291.       printf("Cannot open %s\n","laplace.mat");
  292.       exit(1);
  293.    }
  294.    
  295.    /* b must be transformed to a matrix */
  296.    
  297.    B = m_get(N,N);
  298.    VEC2MAT(out,B);
  299.    m_save(fp,B,"sol");  /* sol is an internal name in MATLAB */
  300.  
  301. }
  302.  
  303.  
  304. void main()
  305. {
  306.    int i;
  307.  
  308.    input("Choose the problem (1=Runge-Kutta, 2=least squares,3=laplace): ",
  309.      "%d",&i);
  310.    switch (i) {
  311.     case 1: tutor_rk4(); break;
  312.     case 2: tutor_ls(); break;
  313.     case 3: tut_lap(); break;
  314.     default: 
  315.       printf(" Wrong value of i (only 1, 2 or 3)\n\n");
  316.       break;
  317.    }
  318.  
  319. }
  320.  
  321.