Back to index

texmacs  1.0.7.15
equations.cpp
Go to the documentation of this file.
00001 
00002 /******************************************************************************
00003 * MODULE     : equations.cpp
00004 * DESCRIPTION: Some tools for solving systems of equations
00005 * COPYRIGHT  : (C) 2004  Henri Lesourd
00006 *******************************************************************************
00007 * This software falls under the GNU general public license version 3 or later.
00008 * It comes WITHOUT ANY WARRANTY WHATSOEVER. For details, see the file LICENSE
00009 * in the root directory or <http://www.gnu.org/licenses/gpl-3.0.html>.
00010 ******************************************************************************/
00011 
00012 #include "point.hpp"
00013 #include "math_util.hpp"
00014 
00015 // Tridiagonal & quasi tridiagonal systems
00016 void
00017 tridiag_solve (array<double> a, array<double> b, array<double> c,
00018               array<point> x, array<point> y, int n)
00019 {
00020   array<double> u(n);
00021   int i;
00022   double li;
00023   li= b[0];
00024   ASSERT (b[0] != 0, "failed tridiag_solve (1)");
00025   x[0]= y[0]/li;
00026   for (i=0; i<n-1; i++) {
00027      u[i]= c[i]/li;
00028      li= b[i+1] - a[i+1]*u[i];
00029      ASSERT (li != 0, "failed tridiag_solve (2)");
00030      x[i+1]= (y[i+1] - a[i+1]*x[i]) / li;
00031   }
00032   for (i=n-2; i>=0; i--) {
00033      x[i]= x[i] - u[i]*x[i+1];
00034   }
00035 }
00036 
00037 void
00038 quasitridiag_solve (array<double> a, array<double> b, array<double> c,
00039                   array<double> u, array<double> v,
00040                   array<point> x, array<point> y, int n)
00041 {
00042   int i;
00043   array<point> z(n), up(n);
00044   for (i=0; i<n; i++) up[i]= as_point (u[i]);
00045   tridiag_solve (a, b, c, x, y, n);
00046   tridiag_solve (a, b, c, z, up, n);
00047   point vx;
00048   vx= v[0]*x[0];
00049   for (i=1; i<n; i++) vx= vx + v[i]*x[i];
00050   double vz;
00051   vz= v[0]*z[0][0];
00052   for (i=1; i<n; i++) vz+= v[i]*z[i][0];
00053   vx= vx / (1+vz);
00054   for (i=0; i<n; i++) {
00055      x[i]= x[i] - z[i][0]*vx;
00056   }
00057 }
00058 
00059 void
00060 xtridiag_solve (array<double> a, array<double> b, array<double> c,
00061               double a0, double a1,
00062               array<point> x, array<point> y, int n)
00063 {
00064   array<double> u(n), v(n);
00065   int i;
00066   for (i=0; i<n; i++) u[i]= v[i]= 0;
00067   u[0]= u[n-1]= 1;
00068   v[0]= a0;
00069   v[n-1]= a1;
00070   b[0]-= a0;
00071   b[n-1]-= a1;
00072   quasitridiag_solve (a, b, c, u, v, x, y, n);
00073   b[0]+= a0;
00074   b[n-1]+= a1;
00075 }