home *** CD-ROM | disk | FTP | other *** search
/ ARM Club 1 / ARM_CLUB_CD.iso / contents / apps / clib / progs / meschach / !Meschach / c / zfunc < prev    next >
Encoding:
Text File  |  1994-08-03  |  4.4 KB  |  213 lines
  1. ,"f");
  2.    
  3.    out->ve[0] = x->ve[1];
  4.    out->ve[1] = - x->ve[0];
  5.    
  6.    return out;
  7. }
  8.  
  9.  
  10. void tutor_rk4()
  11. {
  12.    VEC        *x;
  13.    VEC        *f();
  14.    double     h, t, t_fin;
  15.    double     rk4();
  16.    
  17.    input("Input initial time: ","%lf",&t);
  18.    input("Input final time: ",  "%lf",&t_fin);
  19.    x = v_get(2);    /* this is the size needed by f() */
  20.    prompter("Input initial state:\n");    x = v_input(VNULL);
  21.    input("Input step size: ",   "%lf",&h);
  22.    
  23.    printf("# At time %g, the state is\n",t);
  24.    v_output(x);
  25.    while (t < t_fin)
  26.    {
  27.       /* you can use t = rk4_var(f,t,x,min(h,t_fin-t)); */
  28.       t = rk4(f,t,x,min(h,t_fin-t));   /* new t is returned */
  29.       printf("# At time %g, the state is\n",t);
  30.       v_output(x);
  31.    }
  32. }
  33.  
  34.  
  35.  
  36.  
  37. #include "matrix2.h"
  38.  
  39. void tutor_ls()
  40. {
  41.    MAT *A, *QR;
  42.    VEC *b, *x, *diag;
  43.    
  44.    /* read in A matrix */
  45.    printf("Input A matrix:\n");
  46.    
  47.    A = m_input(MNULL);     /* A has whatever size is input */
  48.    
  49.    if ( A->m < A->n )
  50.    {
  51.       printf("Need m >= n to obtain least squares fit\n");
  52.       exit(0);
  53.    }
  54.    printf("# A =\n");       m_output(A);
  55.    diag = v_get(A->m);
  56.    /* QR is to be the QR factorisation of A */
  57.    QR = m_copy(A,MNULL);
  58.    QRfactor(QR,diag);   
  59.    /* read in b vector */
  60.    printf("Input b vector:\n");
  61.    b = v_get(A->m);
  62.    b = v_input(b);
  63.    printf("# b =\n");       v_output(b);
  64.    
  65.    /* solve for x */
  66.    x = QRsolve(QR,diag,b,VNULL);
  67.    printf("Vector of best fit parameters is\n");
  68.    v_output(x);
  69.    /* ... and work out norm of errors... */
  70.    printf("||A*x-b|| = %g\n",
  71.       v_norm2(v_sub(mv_mlt(A,x,VNULL),b,VNULL)));
  72. }
  73.  
  74.  
  75. #include "iter.h"
  76.  
  77.  
  78. #define N 50
  79. #define VEC2MAT(v,m)  vm_move((v),0,(m),0,0,N,N);
  80.  
  81. #define PI 3.141592653589793116
  82. #define index(i,j) (N*((i)-1)+(j)-1)
  83.  
  84. /* right hand side function (for generating b) */
  85. double f1(x,y)
  86. double x,y;
  87. {
  88.   /* return 2.0*PI*PI*sin(PI*x)*sin(PI*y); */
  89.    return exp(x*y);
  90. }
  91.  
  92. /* discrete laplacian */
  93. SPMAT *laplacian(A)
  94. SPMAT *A;
  95. {
  96.    Real h;
  97.    int i,j;
  98.    
  99.    if (!A)
  100.      A = sp_get(N*N,N*N,5);
  101.  
  102.    for ( i = 1; i <= N; i++ )
  103.      for ( j = 1; j <= N; j++ )
  104.      {
  105.         if ( i < N )
  106.       sp_set_val(A,index(i,j),index(i+1,j),-1.0);
  107.         if ( i > 1 )
  108.       sp_set_val(A,index(i,j),index(i-1,j),-1.0);
  109.         if ( j < N )
  110.       sp_set_val(A,index(i,j),index(i,j+1),-1.0);
  111.         if ( j > 1 )
  112.       sp_set_val(A,index(i,j),index(i,j-1),-1.0);
  113.         sp_set_val(A,index(i,j),index(i,j),4.0);
  114.      }
  115.    return A;
  116. }
  117.  
  118. /* generating right hand side */
  119. VEC *rhs_lap(b)
  120. VEC *b;
  121. {
  122.    Real h,h2,x,y;
  123.    int i,j;
  124.    
  125.    if (!b)
  126.      b = v_get(N*N);
  127.  
  128.    h = 1.0/(N+1);      /* for a unit square */
  129.    h2 = h*h;
  130.    x = 0.0;
  131.    for ( i = 1; i <= N; i++ ) {
  132.       x += h;
  133.       y = 0.0;
  134.      for ( j = 1; j <= N; j++ ) {
  135.     y += h;
  136.     b->ve[index(i,j)] = h2*f1(x,y);
  137.      }
  138.    }
  139.    return b;
  140. }
  141.    
  142. void tut_lap()
  143. {
  144.    SPMAT *A, *LLT;
  145.    VEC *b, *out, *x;
  146.    MAT *B;
  147.    int num_steps;
  148.    FILE *fp;
  149.  
  150.    A = sp_get(N*N,N*N,5);
  151.    b = v_get(N*N);
  152.  
  153.    laplacian(A);
  154.    LLT = sp_copy(A);
  155.    spICHfactor(LLT);
  156.  
  157.    out = v_get(A->m);
  158.    x = v_get(A->m);
  159.  
  160.    rhs_lap(b);   /* new rhs */
  161.    iter_spcg(A,LLT,b,1e-6,out,1000,&num_steps);
  162.    printf("Number of iterations = %d\n",num_steps);
  163.  
  164.    /* save b as a MATLAB matrix */
  165.  
  166.    fp = fopen("laplace.mat","w");  /* b will be saved in laplace.mat */
  167.    if (fp == NULL) {
  168.       printf("Cannot open %s\n","laplace.mat");
  169.       exit(1);
  170.    }
  171.    
  172.    /* b must be transformed to a matrix */
  173.    
  174.    B = m_get(N,N);
  175.    VEC2MAT(out,B);
  176.    m_save(fp,B,"sol");  /* sol is an internal name in MATLAB */
  177.  
  178. }
  179.  
  180.  
  181. void main()
  182. {
  183.    int i;
  184.  
  185.    input("Choose the problem (1=Runge-Kutta, 2=least squares,3=laplace): ",
  186.      "%d",&i);
  187.    switch (i) {
  188.     case 1: tutor_rk4(); break;
  189.     case 2: tutor_ls(); break;
  190.     case 3: tut_lap(); break;
  191.     default: 
  192.       printf(" Wrong value of i (only 1, 2 or 3)\n\n");
  193.       break;
  194.    }
  195.  
  196. }
  197.  
  198. FileDataŵupdateq
  199. EÿÿÿÈ(?Jd
  200. /**************************************************************************
  201. **
  202. ** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  203. **
  204. **                 Meschach Library
  205. ** 
  206. ** This Meschach Library is provided "as is" without any express 
  207. ** or implied warranty of any kind with respect to this software. 
  208. ** In particular the authors shall not be liable for any direct, 
  209. ** indirect, special, incidental or consequential damages arising 
  210. ** in any way from use of the software.
  211. ** 
  212. ** Everyone is granted permission to copy, modify and redistribute this
  213. ** Meschach Library, provided:
  214. **  1.  All copies con