ex5.c

multipleshooting to find the roots of the non-linear functions

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include "math_util.h"
#define M 10.0
#define D 5.0
#define L 5.0
#define In 12.0
long double alpha;
main()
{
   char file_name[12];
   FILE *output_file, *input_file;
   extern long double alpha;
   int itt, n, m, ms, i, j, k, ierr;
   long double tend, **yout, **zout, *time, tol;

   void dyn( long double *, long double *, long double *, long double * );
   void dynbc( long double *, long double *, long double * );

   void dynjac( long double *, long double *, long double **, 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 = 12;
   m = 4;
   ms = 21;
   time = a1d_allo_dbl( ms );
   yout = a2d_allo_dbl( ms, n );
   zout = a2d_allo_dbl( ms, m );
   tol = 1.0e-3;
   tend = 5.0;
   alpha =1.0e-3;
     
   mat_null( yout, ms, n );
   mat_null( zout, ms, m );
   for( i = 0; i < ms; i++ ) {
      time[i] = i*tend/(ms-1);
      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]);   
   }
   
   itt=0;
   while( itt < 9) {
      sprintf(file_name,"tpbvl_%d",itt++);
      output_file = fopen(file_name,"w");
      if( output_file == NULL ) {
         printf("ERROR unable to open output file\n");
         exit(1);
      }
      ierr = mshoot_dae( yout, zout, time, n, m, ms, tol, dyn, dynjac, dynbc, dyndbc );
      if( ierr != 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");
      }
      i=ms-1;
      printf("%Le %Le %Le %Le %Le %Le\n",time[i],yout[i][0],yout[i][2],yout[i][4],zout[i][0],zout[i][1]);   
    
      fclose( output_file );
   
      alpha *= 5.0;   
   }
}
void dynbc( sy0, syms, r )
long double *sy0, *syms, *r;
{
   extern long double alpha;
   
   r[0] = sy0[0];
   r[1] = sy0[1];
   r[2] = sy0[2];
   r[3] = sy0[3];
   r[4] = sy0[4];
   r[5] = sy0[5];
   r[6] = syms[6]-alpha*(syms[0]-4.0);
   r[7] = syms[7]-alpha*syms[1];
   r[8] = syms[8]-alpha*(syms[2]-4.0);
   r[9] = syms[9]-alpha*syms[3];
   r[10] = syms[10]-alpha*(syms[4]-0.785398164);
   r[11] = syms[11]-alpha*syms[5];
}
void dyn( y, z, f, g )
long double *y, *z, *f, *g;
{
   long double s4, c4;
   
   s4 = sin(y[4]);
   c4 = cos(y[4]);
   
   f[0] =  y[1];
   f[1] =  ((z[0]+z[2])*c4-(z[1]+z[3])*s4)/M;
   f[2] =  y[3];
   f[3] =  ((z[0]+z[2])*s4+(z[1]+z[3])*c4)/M;
   f[4] =  y[5];
   f[5] =  ((z[0]+z[2])*D+(-z[1]+z[3])*L)/In;
   f[6] =  0.0;
   f[7] = -y[6];
   f[8] =  0.0;
   f[9] = -y[8];
   f[10] =  y[7]*f[3]-y[9]*f[1];
   f[11] = -y[10];
   g[0] =  y[7]*c4/M + y[9]*s4/M + y[11]*D/In + z[0];
   g[1] = -y[7]*s4/M + y[9]*c4/M - y[11]*L/In + z[1];
   g[2] =  y[7]*c4/M + y[9]*s4/M + y[11]*D/In + z[2];
   g[3] = -y[7]*s4/M + y[9]*c4/M + y[11]*L/In + z[3];
}
void dyndbc( sy0, syms, a, b )
long double *sy0, *syms, **a, **b;
{
   extern long double alpha;
   int i;
   
   mat_null( a, 12, 12 );
   mat_null( b, 12, 12 );
   
   for(i=0;i<6;i++) {
      a[i][i]=1.0;
      b[i+6][i]=-alpha;
      b[i+6][i+6]=1.0;
   }
}void dynjac( y, z, fy, fz, gy, gz )
long double *y, *z, **fy, **fz, **gy, **gz;
{
   long double s4, c4;
   int n = 12, m = 4;
   void mat_null( long double **, int, int );
     
   mat_null( fy, n, n );
   mat_null( fz, n, m );
   mat_null( gy, m, n );
   mat_null( gz, m, m );   
   
   s4 = sin(y[4]);
   c4 = cos(y[4]);
   
   fy[0][1] =  1.0;
   fy[1][4] = -((z[0]+z[2])*s4+(z[1]+z[3])*c4)/M;
   fy[2][3] =  1.0;
   fy[3][4] =  ((z[0]+z[2])*c4-(z[1]+z[3])*s4)/M;
   fy[4][5] =  1.0;
   fy[7][6] = -1.0;
   fy[9][8] = -1.0;
   fy[10][4] =  y[7]*(((z[0]+z[2])*c4-(z[1]+z[3])*s4)/M)+y[9]*(((z[0]+z[2])*s4+(z[1]+z[3])*c4)/M);
   fy[10][7] =  ((z[0]+z[2])*s4+(z[1]+z[3])*c4)/M;
   fy[10][9] = -((z[0]+z[2])*c4-(z[1]+z[3])*s4)/M;
   fy[11][10] = -1.0;
   
   fz[1][0] = c4/M;fz[1][2] = fz[1][0];
   fz[1][1] = -s4/M;fz[1][3] = fz[1][1];
   fz[3][0] = s4/M;fz[3][2] = fz[3][0];
   fz[3][1] = fz[1][0];fz[3][3] = fz[1][0];
   fz[5][0] = D/In;fz[5][2] = fz[5][0];
   fz[5][1] = -L/In;
   fz[5][3] = L/In;
   fz[10][0] = y[7]*s4/M-y[9]*c4/M;fz[10][2]=fz[10][0];
   fz[10][1] = y[7]*c4/M+y[9]*s4/M;fz[10][3]=fz[10][1];
   
   gy[0][4]=-y[7]*s4/M + y[9]*c4/M;
   gy[0][7]=c4/M;
   gy[0][9]=s4/M;
   gy[0][11]=D/In;
   gy[1][4]=-y[7]*c4/M - y[9]*s4/M;
   gy[1][7]=-s4/M;
   gy[1][9]=c4/M;
   gy[1][11]=-L/In;
   gy[2][4]=-y[7]*s4/M + y[9]*c4/M;
   gy[2][7]=c4/M;
   gy[2][9]=s4/M;
   gy[2][11]=D/In;
   gy[3][4]=-y[7]*c4/M - y[9]*s4/M;
   gy[3][7]=-s4/M;
   gy[3][9]=c4/M;
   gy[3][11]=L/In;
   
   gz[0][0]=1.0;gz[1][1]=1.0;gz[2][2]=1.0;gz[3][3]=1.0;
   
}