ex6.c

multipleshooting to find the roots of the non-linear functions

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include "math_util.h"
/*

  This problem is from Bryson, A. and Ho, Y., Applied Optimal Control, pages 52-55
  Minimum drag nose shape in hypersonic flow
  
*/
main()
{
   FILE *output_file;
   int n, m, ms, i, j, ierr;
   long double tend, y0, y1, y2, y3, z0, **yout, **zout, *time, tol, T, dt, x;

   void dyn( long double *, long double *, long double *, long double * );
   void dynbc( long double *, long double *, long double * );
   int mshoot_dae( long double **, long double **, long double *, int, int, int, long double, 
                   void(*)(), void(*)(), void(*)(), void(*)() );
   
   void dynjac( long double *, long double *, long double **, long double **, long double **, long double ** );
   void dyndbc( long double *, long double *, long double **, long double ** );

   n = 2;
   m = 1;
   ms = 101;
   time = a1d_allo_dbl( ms );
   yout = a2d_allo_dbl( ms, n );
   zout = a2d_allo_dbl( ms, m );
   tol = 1.0e-4;
   tend = 4.0;
   
   output_file = fopen("ex6.dat","w");
   if( output_file == NULL ) {
      printf("ERROR unable to open output file\n");
      exit(1);
   }
/*

   Setup initial guess

*/   
   y0 = 1;
   y1 = 0.25;
   z0 = 0.3;
  
   T = tend;
   dt = 1.0/(ms-1);
   for( i = 0; i < ms; i++ ) {
      x = pow(10.0,i*dt);
      time[ms-1-i]=T+T*(1.0-x)/9.0;
   }

   for( i = 0; i < ms; i++ ) {
      yout[i][0] = y0;
      yout[i][1] = y1;
      zout[i][0] = z0;
      printf("%Le %Le %Le %Le\n",time[i],yout[i][0],yout[i][1],zout[i][0]);
   }
    
   ierr = mshoot_dae( yout, zout, time, n, m, ms, tol, dyn, dynjac, dynbc, dyndbc );

   x = zout[0][0];
   dt = 3.0+10.0*x*x+17.0*x*x*x*x+2.0*x*x*x*x*x*x+4.0*x*x*x*x*log(1.0/x);
   dt *= (x*x)/((1.0+x*x)*(1.0+x*x)*(1.0+x*x)*(1.0+x*x));
   
   printf("Cd = %Le\n",dt); /* Coefficient of drag */
   
   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]);
      dt = 2.0*yout[i][0]*zout[i][0];
      dt *= zout[i][0]*zout[i][0];
      dt /= ((1.0+zout[i][0]*zout[i][0])*(1.0+zout[i][0]*zout[i][0]));
      dt *= -1.0;
      fprintf(output_file,"%Le\t",dt);
      fprintf(output_file,"\n");
   }
    
   fclose( output_file );
   exit(0);
}
void dynbc( sy0, syms, r )
long double *sy0, *syms, *r;
{
   r[0] = sy0[0]-1.0;
   r[1] = syms[1]-syms[0];
}
void dyndbc( sy0, syms, a, b )
long double *sy0, *syms, **a, **b;
{
   void mat_null( long double **, int, int );
   
   mat_null( a, 2, 2 );
   mat_null( b, 2, 2 );
   
   a[0][0] = 1.0;
   b[1][0] = -1.0;
   b[1][1] = 1.0;
}
void dyn( y, z, f, g )
long double *y, *z, *f, *g;
{
   long double u;
   
   u = z[0];
   f[0] = -u;
   f[1] = -(u*u*u)/(1.0+u*u);
   g[0] = (y[0]*u*u*(3.0+u*u))/((1.0+u*u)*(1.0+u*u)) - y[1];
}


void dynjac( y, z, fy, fz, gy, gz )
long double *y, *z, **fy, **fz, **gy, **gz;
{
   void mat_null( long double **, int, int );
   int n = 2, m = 1;
   long double u, d;
   
   u = z[0];
   d = 1.0+u*u;
   
   mat_null( fy, n, n );
   mat_null( fz, n, m );
   mat_null( gy, m, n );
   mat_null( gz, m, m );
   
   
   fz[0][0] = -1.0;
   fz[1][0] = -3.0*u*u/(1.0+u*u) + (2.0*u*u*u*u)/(d*d);
   
   gy[0][0] = (u*u*(3.0+u*u))/(d*d);
   gy[0][1] = -1.0;
   
   gz[0][0] = (2.0*y[0]*u*(3.0+u*u))/(d*d) + (2.0*y[0]*u*u*u)/(d*d) - (4.0*y[0]*u*u*u*(3.0+u*u))/(d*d*d);
}