ex9a.c

multipleshooting to find the roots of the non-linear functions

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include "math_util.h"
long double alpha;
main()
{
   FILE *output_file, *input_file;
   extern long double alpha;
   int n, m, ms, i, j, k, ierr;
   long double tend, y0, y1, y2, y3, z0, z1, z2, **yout, **zout, *time, tol;

   void dyn( long double *, long double *, long double *, long double * );
   void dynjac( long double *, long double *, long double **, long double **, long double **, long double ** );
   void dynbc( long double *, long double *, long double * );
   void dyndbc( long double *, long double *, long double **, long double ** );
   int mshoot_dae( long double **, long double **, long double *, int, int, int, long double, 
                   void(*)(), void(*)(), void(*)(), void(*)() );
   
   n = 4;
   m = 3;
   ms = 51;
   time = a1d_allo_dbl( ms );
   yout = a2d_allo_dbl( ms, n );
   zout = a2d_allo_dbl( ms, m );
   tol = 1.0e-4;
   tend = 10.0;
   alpha = 0.1;
   
   output_file = fopen("ex9a.dat","w");
   if( output_file == NULL ) {
      printf("ERROR unable to open output file\n");
      exit(1);
   }
   y0 = 0.0;
   y1 = 0.0;
   y2 = 0.1;
   y3 = -0.025;
   z0 = 1.0;
   z1 = 1.0;
   z2 = 0.0;

   for( i = 0; i < ms; i++ ) {
      time[i] = i*tend/(ms-1);
      yout[i][0] = 3.66-(4.66/tend)*time[i];
      yout[i][1] = -1.86+(2.86/tend)*time[i];
      yout[i][2] = y2;
      yout[i][3] = y3;
      zout[i][0] = 0.5;zout[i][1] = sqrt(1.0-zout[i][0]*zout[i][0]);
      zout[i][2] = 1.0;
      printf("%Le %Le %Le %Le %Le %Le\n",time[i],yout[i][0],yout[i][1],yout[i][2],yout[i][3],zout[i][0]);
   }
   for( k = 0; k < 11; k++ ) {
       ierr = mshoot_dae( yout, zout, time, n, m, ms, tol, dyn, dynjac, dynbc, dyndbc );
       alpha += 0.1;
       if( alpha >=1.0 )
          break;
   }
   for( i = 0; i < ms; i++ ) {
      fprintf(output_file,"%Le\t",time[i]);
      for( j = 0; j < n; j++ )
         fprintf(output_file,"%Le\t",yout[i][j]);
      for( j = 0; j < m; j++ )
         fprintf(output_file,"%Le\t",zout[i][j]);
      fprintf(output_file,"\n");
   }
   pause();
    
   fclose( output_file );    
}
void dynbc( sy0, syms, r )
long double *sy0, *syms, *r;
{
   extern long double alpha;
   r[0] = sy0[0]-3.66;
   r[1] = sy0[1]+1.86;
   r[2] = syms[2]-syms[0];
   r[3] = syms[3]-syms[1];
}
void dyndbc( sy0, syms, a, b )
long double *sy0, *syms, **a, **b;
{
   extern long double alpha;
   mat_null( a, 4, 4 );
   mat_null( b, 4, 4 );
   a[0][0] = 1.0;
   a[1][1] = 1.0;
   b[2][2] = 1.0;b[2][0] = -1.0;
   b[3][3] = 1.0;b[3][1] = -1.0;
}
void dyn( y, z, f, g )
long double *y, *z, *f, *g;
{
   extern long double alpha;
   
   f[0] = -alpha*y[1]+z[0];
   f[1] = z[1];
   f[2] = 0.0;
   f[3] = alpha*y[2];
   g[0] = y[2]+2.0*z[0]*z[2];
   g[1] = y[3]+2.0*z[1]*z[2];
   g[2] = z[0]*z[0]+z[1]*z[1]-1.0;
}
void dynjac( y, z, fy, fz, gy, gz )
long double *y, *z, **fy, **fz, **gy, **gz;
{
   extern long double alpha;
   int n = 4, m = 3;
   
   mat_null( fy, n, n );
   mat_null( fz, n, m );
   mat_null( gy, m, n );
   mat_null( gz, m, m );
   
   fy[0][1] = -alpha;
   fy[3][2] = alpha;
   
   fz[0][0] = 1.0;fz[1][1] = 1.0;
   
   gy[0][2] = 1.0;gy[1][3] = 1.0;
   
   gz[0][0] = 2.0*z[2];
   gz[0][2] = 2.0*z[0];
   gz[1][1] = 2.0*z[2];
   gz[1][2] = 2.0*z[1];
   gz[2][0] = 2.0*z[0];
   gz[2][1] = 2.0*z[1];
}