/* v2a.c value function integration 5 November, 1999 ver 9.6 C. Rowat cir20@cam.ac.uk cc -o v2a.exe v2a.c -I/opt/nagc/include -R/opt/nagc -L/opt/nagc -lnagc -lm CUS CPU limit commands: ulimit -a displays, ulimit -St XX sets, howmuch displays usage SIGNALLING PHASE TRANSITIONS The global variable sflag is used to track whether players are in the interior or on the corner; this is managed by sflagger(). sflag is set before the while loop in normalpath() and then updated in the gbothint() and g1int() routines. Previously sflag had only been set immediately inside the while loop; this allowed sflag to be used inside g1int() to identify the cornered player (without the global variable now used). The problem with that approach is that the d02 routines calculated to +/- TOL on z. When d02 returned a corner, for example, sflagger() would not necessarily identify it as such as the v[i] value might still indicate interior: sflagger() was more sensitive than the d02 routines. The use of (> and <) +TOL elements in sflagger() as well as the conventions (>= and <) and (> and <) was also attempted. SWITCHES -b: uses bdf (implicit) NAG routine d02ejc when d02cjc -> NE_TOL_PROGRESS bdf=0: no bdf used (default); bdf=1: bdf if non-bdf fails (switch); bdf=2: non-bdf has failed (temporary); Attempts to continue with d02ejc from NE_TOL_PROGRESS error of d02cjc failed -g: graphic generates LaTeX picture from gridsearch() -l: calculates downward sloping linear path parameters -p: generates .path logfile -r: reverses direction of integration [incomplete implementation] -v: verbose log; automatically implies -p FLAG CODES eflag sflag uflag pcflag 0 gridsearch well cond'd 1 [both explode] both interior - poorly cond'd 2 [0 explodes] 0 interior, 1 corner singshoot 3 [1 explodes] 0 corner, 1 interior singtrap 4 non-invertibility hit both corner specpath 5 [limit of integration hit: fail.code deals] section 6 corner survival 7 [both collapse] 8 [0 collapses] 9 [1 collapses] 10 o.d.e denominator becomes small: advance warning of non-invertibility? 11 quasi non-invertible 12 V > 0: player 0 13 V > 0: player 1 14 [L-infinity norm exceeds TCOND: poor conditioning] 15 V > 0: both players 16 [as 10 but also badly conditioned] 17 [as 12 but also badly conditioned] 18 [as 13 but also badly conditioned] 19 [as 15 but also badly conditioned] 20 non-invertibility and |ZDOT| < TOL [] in the eflags indicate no longer active eflags. A BRIEF HISTORY OF VERSION INNOVATIONS v2a: poorly conditioned paths: treated as normal but pcflag'ed if path EVER is poorly conditioned. This done to avoid the problem of reaching P-C before non-inv'y nomenclature change: v(alue function) 2(player) a(symmetric): asymmetric players v13: knownlinear() cut: superseded by specpath() generalising code to deal with asymmetric players have generalised initial conditions to include any zstart (not nec'y =0.) section() allows sections (1-d grids) to be used as starting points. specpath() allows specific path to be integrated along by spec'g ul and ll as v[] quasi non-invert checks moved into gbothint() and g1int(). Allow poorly conditioned? repeated opening and closing of output files replaced with fflush() commands. allows singtrap() to step either symmetrically or fixing 1 player's init conds traps individual points on the non-invertibility - V > 0 border traces singularity locus v12: allows internal steps to report to be controlled from datafile as well splits .log output into: .path : contains (z,v0,v1) path information .result: contains conclusions on paths reads starting grid information from datafile to allow more flexibility does less in main() looks for more "non-grid" solutions through the singularity allows backwards integration to recover initial conditions; achieved by soft fail handling at zend v10: repairs sflag==4 routine (1/2/99) v9: condenses code by adding new functions v8: updated checking; better integration of linear and general runs v7: reads parameters from data file; different loop structure v3: FEW refinements used (v. those implicitly based on symmetric case) I have avoided use of pow(double,double) as it is slow Trying to prevent first call to a g() returning -1; trap before entry This code uses the following NAG C Mark 5 routines: a00aac: NAG implementation details c02agc: modified Laguerre's method polynomial root finder c05adc: root finder for continuous function of one variable c05tbc: hybrid Powell method root finder for nonlinear system d02cjc: Adam's method IVP ode integration */ #include #include #include #include #include #include #include #include #include #include #include /*A, B, S_B, B_HAT, S_B_HAT only useful in symmetric case*/ #define A (2.*nu[0]*(beta*zeta[0]-xi[0]-xi[1])/(3.*nu[0]+beta*(beta+delta[0]))) #define B ((3.*nu[0]*zeta[0]+xi[0]+xi[1]*(beta+delta[0]))/(3.*nu[0]+beta*(beta+delta[0]))) #define S_B ((2.*beta+delta[0]-sqrt((2.*beta+delta[0])*(2.*beta+delta[0])+12.*nu[0]))/3.) #define ZDOT (((v[0]>-2.*xi[0])?(xi[0]+v[0]/2.):0)+((v[1]>-2.*xi[1])?(xi[1]+v[1]/2.):0)-beta*z) #define VDENOM (ZDOT*ZDOT-(v[0]*v[1])/4.) #define TOL 1E-8 /*machine zero ~ 1E-16=X02AJC; NAG: use sqrt() for TOL*/ #define TCOND 1E10 /*condition number: Judd text recommendation*/ #define B_HAT ((nu[0]*zeta[0]+xi[0]*(beta+delta[0]))/(nu[0]+beta*(beta+delta[0]))) #define S_B_HAT ((2.*beta+delta[0]-sqrt((2.*beta+delta[0])*(2.*beta+delta[0])+4.*nu[0]))/2.) /*C and Q macros used for v^dot = r + lambda * q calculations along singularity locus*/ #define C1 (sqrt(v1/v0)*(v0<0?1.:-1.)) #define C2 1. #define Q1 (sqrt(v0/v1)*(v0<0?1.:-1.)) #define Q2 1. static void gridsearch(double zstart); static void singshoot(); static void singtrap(); static void specpath(double zstart); static void section(double zstart); static void singloc(Integer neq, double args[], double fout[], Integer *userflag, Nag_User *comm); static int typeset(); static void calc_slopes(double z, double v0, double v1, double *lam_ptr); static void calc_devns(double z, double v0, double v1, double vp[], double *dz_ptr); static void normalpath(double z, double v[], double *state_ptr); static double setpara(); static void out(Integer neq, double *zsol, double v[], Nag_User *comm); static double gbothint(Integer neq, double z, double v[], Nag_User *comm); static double g1int(Integer neq, double z, double v[], Nag_User *comm); static int gevent(Integer code); static void vprime(Integer neq, double z, double v[], double vp[], Nag_User *comm); static double deriv(int i,double z, double v[],int sysflag); static double kroot(double z); static void sflagger(double z, double v[]); static int lpfiler(Integer i, Integer j, Integer lastdot); static void jacob1(Integer neq, double z, double v[], double pw[], Nag_User *comm); static void jacobboth(Integer neq, double z, double v[], double pw[], Nag_User *comm); static void lpsetup(); static void calclinear(); static void lincalc(Integer neq, double args[], double fout[], Integer *userflag, Nag_User *comm); static void lincalc2(Integer neq, double args[], double fout[], Integer *userflag, Nag_User *comm); struct user { double zend, h; Integer k; }; double beta, delta[2], nu[2], xi[2], zeta[2]; /*game parameters*/ double zend, ul, ll, ul2, ll2; /*path parameters*/ char scen_name[40]; /*scenario name*/ double ki; int whichi; /*passes player info from normalpath() to kroot*/ int sflag; /*which equations of motion in use*/ int whocorn; /*"who has cornered" flag*/ int eflag; /*event flag*/ int pcflag; /*poor conditioning flag*/ Integer dsign; /*tracks denom sign for sing'y*/ int uflag,initcond, int_step; /*programme option variables*/ FILE *logr, *logp, *loglp; /*.result, .path and .lp logfile pointers*/ int verbose=0, lp_on=0, symmetric=0, reverse=0, bdf=0, path=0, linear=0;/*switches*/ main(int argc, char *argv[]) /********************************************************************************\ * Reads parameters before directing programme to: * * gridsearch(): walk over a grid of initial conditions; * * singshoot(): trace the singularity locus, integrating off of it; * * singtrap(): trap a particular point on the singularity. * \********************************************************************************/ { double zstart; /*initial z for integration*/ time_t t=time(NULL); /*to count seconds*/ int i; /*parameter display*/ a00aac(); printf(" machine zero is %20.18f\n",X02AJC); logr = fopen("v2a.result","a"); fprintf(logr,"TWO ASYMMETRIC PLAYER VALUE FUNCTION NUMERICAL INTEGRATION\n"); fprintf(logr,"time: %ld\n\n",t); fflush(logr); zstart=setpara(); /*sets parameters*/ fprintf(logr,"Scenario: %s\n",scen_name); fprintf(logr,"SYSTEM PARAMETER\n beta=%13.8f, zend=%8.5f, uflag=%d\n",beta,zend,uflag); fprintf(logr,"PLAYER PARAMETERS\n"); for(i=0;i<2;++i) { fprintf(logr,"Player %d: delta=%13.8f, nu=%13.8f, xi=%13.8f, zeta=%13.8f\n", i,delta[i],nu[i],xi[i],zeta[i]); } fprintf(logr,"\n");fflush(logr); while ((argc>1)&&(argv[1][0]=='-')) /*parses command line*/ { switch (argv[1][1]) { case 'b': /*retries path with BDF routine d02ejc if d02cjc -> NAG failure*/ bdf=1; break; case 'g': /*generates crude .lp output file*/ lp_on=1; break; case 'l': linear=1; break; case 'p': path=1; break; case 'r': reverse=1; break; case 'v': /*verbose logging*/ verbose=1; path=1; break; } ++argv; --argc; } if (path==1) logp=fopen("v2a.path","a"); if((lp_on==1)&&(uflag==0)) /*only .lp for gridsearch*/ { loglp=fopen("v2a.lp","w"); /*replacing .lp file rather than appending*/ } if (linear==1) calclinear(); if (uflag==0) gridsearch(zstart); else if (uflag==2) singshoot(); else if (uflag==3) singtrap(); else if (uflag==4) specpath(zstart); else if (uflag==5) section(zstart); fclose(logr); if(loglp!=NULL) fclose(loglp); if(logp!=NULL) fclose(logp); exit(EXIT_SUCCESS); } static void gridsearch(double zstart) /********************************************************************************\ * The basic function for generating the phase plot: walks over a grid of initial * * conditions. * \********************************************************************************/ { double v[2], last_state[3]; Integer i,j,h,lastdot=-1; if(verbose==1) {fprintf(logp,"gridsearch() entered\n"); fflush(logp);} if(lp_on==1) lpsetup(); /*writes preamble to .lp file*/ for (i=0;i<=initcond;++i) { if(symmetric==1) /*this formatting only if 1/2 paths tracked*/ { for (h=1;h<=i;++h) fprintf(logr,", "); /*prints the leading commas*/ fflush(logr); } /*lastdot=-1; here -> dots on diagonal line*/ for (j=(symmetric==1?i:1);j<=initcond;++j) /*only half if symmetry*/ { v[0]=ll+(ul-ll)*((double)i/(double)initcond); v[1]=ll+(ul-ll)*((double)j/(double)initcond); normalpath(zstart,v,last_state); if (bdf==2) { if(verbose==1) {fprintf(logp,"recalling normalpath() with bdf==2\n"); fflush(logp);} /*array names automatically pointers; therefore inherit v[] from normalpath()*/ v[0]=ll+(ul-ll)*((double)i/(double)initcond); v[1]=ll+(ul-ll)*((double)j/(double)initcond); normalpath(zstart,v,last_state); } if (path==1) {fprintf(logp,"\n"); fflush(logp);} /*outputting results to .lp file*/ if(lp_on==1) lastdot=lpfiler(i,j,lastdot); } /*end of j loop*/ fprintf(logr,"\n"); fflush(logr); } /*end of i loop*/ if(lp_on==1) fprintf(loglp,"\\end{picture}\n");fflush(loglp); } static void singshoot() /********************************************************************************\ * This routine traces the singularity locus by fixing v[0] points (starting at * * -2\xi and continuing either to A if the players are symmetric or to v[1]=-2\xi * * otherwise) and then calculating the v[1] and z values that satisfy the * * equations of the singularity locus. After each (v[0],v[1],z) point on the * * locus is found forward and backward integration then allows the path to be * * calculated. The first steps in this use the slopes calculated by Taylor's * * expansion. * \********************************************************************************/ { /* * VARIABLES: i, j & k index loops; args[0] is used to pass v1 and args[1] to pass * z; vp[] are the slopes, lambda are the elements used in the slope calculation * (v^{dot}=r + lambda*q) and dz is the z dev'n before being able to resume normal * integration. */ Integer i,j,k; double args[2], fout[2], v0, v1, z, v[2], vp[2], last_state[3], lambda[2], dz[2], zstart; double endpoint; /*v0 coordiniate of other end of locus*/ static NagError fail; Nag_User comm; double R1, R2; /*assigned once z has a value; care: also defined locally in calc_slopes*/ if(verbose==1) {fprintf(logp,"singshoot() entered\n"); fflush(logp);} fprintf(logr,"\nPoints along the singular locus\n"); fprintf(logr," v[0] v[1] z fout[0] fout[1]\n"); fflush(logr); /*calculates v0 coordinate of end of singularity locus*/ if(symmetric==1) endpoint=A; else { /*seeds guess; as below but agent scripts reversed. args[0] serving as v[0]*/ args[0]=2.*(nu[0]*xi[1]*(beta*zeta[0]-xi[0])+nu[1]*(beta*zeta[1]-xi[0])) /((beta*(beta+delta[0])+xi[1]*(nu[0]+beta*(beta+delta[1]))+nu[1]+nu[0])*xi[1]+nu[1]); args[1]=(xi[1]*(nu[1]*zeta[1]+nu[0]*zeta[0]*(1.+xi[1]) +(beta+delta[0]+(beta+delta[1])*xi[1])*xi[0])+nu[1]*zeta[1]) /((beta*(beta+delta[0])+xi[1]*(nu[0]+beta*(beta+delta[1]))+nu[1]+nu[0])*xi[1]+nu[1]); comm.p=(Pointer)(-2.*xi[1]); /*initialises v[1] as well*/ c05tbc(2,args,fout,singloc,TOL,&comm,&fail); endpoint=args[0]; } /*seeds initial guesses for v[1] and z, respectively based on sing'y locus eq'ns w/o sqrt()*/ v1=args[0]=2.*(nu[1]*xi[0]*(beta*zeta[1]-xi[1])+nu[0]*(beta*zeta[0]-xi[1])) /((beta*(beta+delta[1])+xi[0]*(nu[1]+beta*(beta+delta[0]))+nu[0]+nu[1])*xi[0]+nu[0]); args[1]=(xi[0]*(nu[0]*zeta[0]+nu[1]*zeta[1]*(1.+xi[0]) +(beta+delta[1]+(beta+delta[0])*xi[0])*xi[1])+nu[0]*zeta[0]) /((beta*(beta+delta[1])+xi[0]*(nu[1]+beta*(beta+delta[0]))+nu[0]+nu[1])*xi[0]+nu[0]); /*steps through the v[0] locus points*/ for(i=0;i<=initcond;++i) { /*step size to traverse from v0=-2xi0 to v1=-2xi1 in initcond steps*/ v0=(double)i/(double)initcond*(endpoint+2.*xi[0])-2.*xi[0]; comm.p=(Pointer)&v0; /*zeros of the nonlinear equations of the singularity locus*/ c05tbc(2,args,fout,singloc,TOL,&comm,&fail); v1=args[0]; /*received after passing*/ fprintf(logr,"%13.5f %13.5f %13.5f %13.5f %13.5f\n",v0,v1,args[1],fout[0],fout[1]); if (fail.code!=NE_NOERROR) fprintf(logr,"%s\n",fail.message); R1=R2=2.*(beta+delta[1]); /*care: assigned in calc_slopes*/ /*gets 2 slopes at sing'y into lambda[] for acceptable points on the singularity locus*/ if ((v1>-2.*xi[1])&&(args[1]>0)) calc_slopes(args[1],v0,v1,lambda); if ((v1>-2.*xi[1])&&(args[1]>0)&&(lambda[0]!=-999.)) /*dismisses complex roots as well*/ { /*steps through both roots*/ for(j=0;j<=1;++j) { /*slope calculations: v_i^{dot} = r_i + lambda * q_i*/ vp[0]=(R1+lambda[j]*Q1)/((v0+v1)/2.+xi[0]+xi[1]-beta*args[1]); vp[1]=(R2+lambda[j]*Q2)/((v0+v1)/2.+xi[0]+xi[1]-beta*args[1]); calc_devns(args[1],v0,v1,vp,dz); /*reads dz[] via pointer*/ for(k=0;k<=1;++k) /*tries both deviations*/ { z=args[1]+dz[k]; v[0]=v0+vp[0]*dz[k]; v[1]=args[0]+vp[1]*dz[k]; while (fabs(VDENOM)<=TOL) /*dz large enough to resume integration?*/ { dz[k]=dz[k]*1.1; z=args[1]+dz[k]; v[0]=v0+vp[0]*dz[k]; v[1]=v1+vp[1]*dz[k]; if(verbose==1) {fprintf(logr,"dz still too small to escape singular region.\n");fflush(logr);} } zstart=setpara(); /*re-initialises zend*/ if (dz[k]>0.) /*forward integration*/ { normalpath(z,v,last_state); if (bdf==2) normalpath(z,v,last_state); fprintf(logr,"\n");fflush(logr); } else /*backward integration*/ { zend=0.; normalpath(z,v,last_state); if (bdf==2) normalpath(z,v,last_state); fprintf(logr,"\n");fflush(logr); } if (path==1) {fprintf(logp,"\n");fflush(logp);} } /*end of k loop*/ } /*end of j loop*/ } /*end of if*/ } /*end of i loop*/ } static void singtrap() /********************************************************************************\ * This function runs a bisection method to find the initial conditions that gen- * * erate a point on the line between a non-invertible path and one violating V>0. * * Triggering this requires setting: uflag=3 * * IT IS NOT CLEAR THAT THIS IS FINDING SINGULARITIES! * \********************************************************************************/ { /* * VARIABLES: The l- refers to lower, the u- to upper. The -start refers to * initial conditions, the -type to the 1,2 or 3 assigned by typeset(). * "input" receives input from the keyboard on "cuttype", whether the bisection * should be along the symmetric axis or with one player fixed. */ static char input[2]; double v[2], lstart, ustart, teststart, last_state[3]; Integer i, ltype, utype, testtype, cuttype; if(verbose==1) {fprintf(logp,"singtrap() entered\n"); fflush(logp);} printf("Type '1' to make a section along the symmetric axis.\n"); printf("Otherwise the section fixes player 0's initial condition.\n"); gets(input); cuttype=atoi(input); ltype=utype=3; /*default does not find the paths of interest*/ v[0]=v[1]=lstart=ll; /*v[0] fixed; only v[1] varies*/ normalpath(0,v,last_state); /*gets path type for lower limit*/ if (bdf==2) normalpath(0,v,last_state); ltype=typeset(); if (path==1) {fprintf(logp,"\n");fflush(logp);} v[1]=ustart=ul; v[0]=((cuttype==1)?ustart:lstart); normalpath(0,v,last_state); /*gets path type for upper limit*/ if (bdf==2) normalpath(0,v,last_state); utype=typeset(); if (path==1) {fprintf(logp,"\n");fflush(logp);} if (ltype==utype) /*no phenomenon*/ { fprintf(logr,"Path type the same at upper and lower limit\n"); fflush(logr); return; } else if (utype*ltype==2) /*phenomenon*/ { while (fabs(lstart-ustart)>2*TOL) { v[1]=teststart=(lstart+ustart)/2.; v[0]=((cuttype==1)?v[1]:ll); normalpath(0,v,last_state); if (bdf==2) normalpath(0,v,last_state); testtype=typeset(); if (path==1) {fprintf(logp,"\n");fflush(logp);} if (testtype==ltype) { lstart=teststart; ltype=testtype; } else if (testtype==utype) { ustart=teststart; utype=testtype; } } } } static void specpath(double zstart) /********************************************************************************\ * Takes ul and ll as v[0] and v[1] initial conditions. Therefore allows specific* * paths to be integrated along. Requires setting uflag=4. * \********************************************************************************/ { double v[2],last_state[3]; if(verbose==1) {fprintf(logp,"specpath() entered\n"); fflush(logp);} v[0]=ll; v[1]=ul; normalpath(zstart,v,last_state); if (bdf==2) { v[0]=ll; v[1]=ul; /*reset as names of arrays are automatically pointers*/ normalpath(zstart,v,last_state); } if((eflag==11)&&(uflag==4)) /*quasi-n-i but want to continue path*/ { if(verbose==1) { fprintf(logp,"Quasi n-i found. Continuing to integrate in opposite direction\n"); fflush(logp); } zend=0.; /*to integrate backward*/ v[0]=last_state[1]; v[1]=last_state[2]; /*seeds call where path left off*/ normalpath(last_state[0],v,last_state); if (bdf==2) normalpath(last_state[0],v,last_state); } } static void section(double zstart) /********************************************************************************\ * Takes ul, ll, ul2 and ll2 as end points for the initial conditions. Then * * integrates along the section with (ul,ul2) and (ll,ll2) as its end points. * * Requires setting uflag=5. * \********************************************************************************/ { double v[2],last_state[3]; Integer i; if(verbose==1) {fprintf(logp,"section() entered\n"); fflush(logp);} for(i=0;i<=initcond;++i) { v[0]=ll+(ul-ll)*((double)i/(double)initcond); v[1]=ll2+(ul2-ll2)*((double)i/(double)initcond); normalpath(zstart,v,last_state); if (bdf==2) normalpath(zstart,v,last_state); if (path==1) {fprintf(logp,"\n"); fflush(logp);} } } static void singloc(Integer neq, double args[], double fout[], Integer *userflag, Nag_User *comm) /********************************************************************************\ * This routine defines the equations whose zeros are to be found. The equations * * are those for non-invertibility and spanning, respectively. Called from sing- * * shoot(). * \********************************************************************************/ { const double G=(*(double *)(comm->p)+args[0])/2.+xi[0]+xi[1]-beta*args[1]; /*dz/dt*/ if(verbose==1) {fprintf(logp,"singloc() entered\n"); fflush(logp);} fout[0]=sqrt(*(double *)(comm->p)*args[0])-2.*((*(double *)(comm->p)+args[0])/2. +xi[0]+xi[1]-beta*args[1]); if(symmetric==1) /*simplified spanning condition for symmetric players*/ fout[1]=G-nu[0]*(args[1]-zeta[0])/(beta+delta[0]); else fout[1]=G*(beta*(2.*G-args[0])+2.*G*delta[0]-args[0]*delta[1]) -args[1]*(2.*G*nu[1]-args[0]*nu[0])+2.*G*nu[1]*zeta[1]-args[0]*nu[0]*zeta[0]; /* if(verbose==1) {fprintf(logp,"v[0]=%13.5f,v[1]=%13.5f,z=%13.5f,fout[0]=%13.5f,fout[1]=%13.5f\n", * *(double *)(comm->p),args[0],args[1],fout[0],fout[1]);fflush(logp);} */ } static int typeset() /********************************************************************************\ * This function returns 1 if an evaluated path hits the truly or the quasi-non- * * invertible locus. * * It returns 2 if the path violates one of the V_i > 0 conditions. * * It returns 3 otherwise. * \********************************************************************************/ { if(verbose==1) {fprintf(logp,"typeset() entered\n"); fflush(logp);} if ((eflag==4)||(eflag==10)||(eflag==11)) return 1; else if ((eflag==12)||(eflag==13)||(eflag==15)) return 2; else return 3; } static void calc_slopes(double z,double v0,double v1,double *lam_ptr) /********************************************************************************\ * Through the singularity there are two slopes, defined by a quadratic equation. * * This routine solves the quadratic and returns the time derivatives. They are * * converted back to z derivatives in singshoot(). We substitute in the deriv- * * before evaluating at the z in question. * \********************************************************************************/ { Complex roos[2]; double coeff[3]; /*coeff[0] for highest order term*/ /* * To substitute into the large quadratic. Note that substitutions have already * occurred on the basis of a_{ij} derivatives, making the coeff[] definitions * somewhat function specific still. The rest of the macros are #defined at the * beginning of this programme. */ const double ZD = ((v0+v1)/2.+xi[0]+xi[1]-beta*z); /*local ZDOT macro*/ const double F11 = ((beta+delta[0])*(ZD+v0/2.)+nu[0]*(z-zeta[0])); const double F12 = ((beta+delta[0])*v0/2.+nu[0]*(z-zeta[0])); const double F13 = (-beta*((beta+delta[0])*v0+2.*nu[0]*(z-zeta[0]))+2.*nu[0]*ZD); const double F21 = ((beta+delta[1])*v1/2.+nu[1]*(z-zeta[1])); const double F22 = ((beta+delta[1])*(ZD+v1/2.)+nu[1]*(z-zeta[1])); const double F23 = (-beta*((beta+delta[1])*v1+2.*nu[1]*(z-zeta[1]))+2.*nu[1]*ZD); const double R1 = 2.*(beta+delta[1]); /*CAREFUL: also assigned in singshoot*/ const double R2 = R1; const double R3 = ZD; coeff[0]=(1.5)*(Q1+Q2); coeff[1]=2.*(R1+R2-beta*R3)+C1*Q2*(R1-F12)+C2*Q1*(R2-F21)-F11-F22; coeff[2]=C1*R1*(R1/2.+R2-beta*R3-F11)+C2*R2*(R1+R2/2.-beta*R3)-C1*(R2*F12+R3*F13) -C2*(R1*F21+R3*F23); if(verbose==1) { fprintf(logp,"calc_slopes() entered\n"); fprintf(logr,"c1=%8.6f,r1=%8.6f,r2=%8.6f,q1=%8.6f,q2=%8.6f,zd=%8.6f\n",C1,R1,R2,Q1,Q2,ZD); fprintf(logr,"f11=%8.6f,f12=%13.11f,f13=%8.6f,f21=%13.11f,f22=%8.6f,f23=%8.6f\n", F11,F12,F13,F21,F22,F23); fprintf(logr,"(2.+C1)*R1=%8.6f\n",(2.+C1)*R1); fprintf(logr,"(2.+Q1)*R2=%8.6f\n",(2.+Q1)*R2); fprintf(logr,"-2.*beta-F11-Q1*F21-C1*F12-F22=%8.6f\n",-2.*beta-F11-Q1*F21-C1*F12-F22); fprintf(logr,"coeff[0]=%8.6f, coeff[1]=%8.6f, coeff[2]=%8.6f\n",coeff[0],coeff[1],coeff[2]); fprintf(logr,"roos[0]=(%8.6f)+(%8.6f)i before c02agc call\n",roos[0].re,roos[0].im); fprintf(logr,"roos[1]=(%8.6f)+(%8.6f)i before c02agc call\n",roos[1].re,roos[1].im); fflush(logr); } c02agc(2,coeff,TRUE,roos,NAGERR_DEFAULT); if(verbose==1) { fprintf(logr,"roos[0]=(%8.6f)+(%8.6f)i after c02agc call\n",roos[0].re,roos[0].im); fprintf(logr,"roos[1]=(%8.6f)+(%8.6f)i after c02agc call\n",roos[1].re,roos[1].im); fflush(logr); } if (!((roos[0].im==0)&&(roos[1].im==0))) { fprintf(logr,"Complex roots detected in slope calculations at v[0]=%13.5f\n",v0); fflush(logr); /*complex roots mean that real paths don't get through*/ *lam_ptr=*(lam_ptr+1)=-999.; /*code for complex roots*/ } else { *lam_ptr=roos[0].re; *(lam_ptr+1)=roos[1].re; /*otherwise return the real roots*/ } return; } static void calc_devns(double z,double v0,double v1,double vp[],double *dz_ptr) /********************************************************************************\ * Before normal NAG integration can begin we need to leave the area in which * * fabs(VDENOM)0.))||(trues[posroot].re<0.)) posroot=i; else if (((trues[i].re>trues[negroot].re)&&(trues[i].re<0.))||(trues[negroot].re>0.)) negroot=i; } /*else complex roots*/ } *dz_ptr=trues[posroot].re; *(dz_ptr+1)=trues[negroot].re; /*what if posroot=negroot=0 still?*/ return; } static void normalpath(double z, double v[], double *state_ptr) /********************************************************************************\ * This routine integrates along a path. It determines the relevant differential * * equations and finally outputs the result to a logfile. * \********************************************************************************/ { Integer neq; double zevent; /*zevent: departure from corner*/ double vp[2]; double k[2]; Nag_User comm; static NagError fail; struct user s; if(verbose==1) {fprintf(logp,"normalpath() entered\n"); fflush(logp);} pcflag=0; /*path initially assumed not poorly conditioned*/ comm.p=(Pointer)&s; s.zend=zend; s.k = int_step; /*internal results to report*/ s.h = (s.zend - z) / (double)(s.k + 1); /*internal result step size*/ neq=2; fail.code=NE_NOERROR; fail.print=FALSE; dsign=0; eflag=0; /*so that subsequent normalpath calls don't inherit the bad eflag*/ /*don't enter while loop for paths starting in corner: they stay there*/ if ((v[0]<-2.*xi[0])&&(v[1]<-2.*xi[1])) eflag=6; sflagger(z,v); /*sets sflag*/ while ((eflag==0)&&(fail.code==NE_NOERROR)) { if(verbose==1) {fprintf(logp,"entering while loop,sflag=%d\n",sflag);fflush(logp);} if (sflag==1) /*both interior*/ { if(verbose==1) {fprintf(logp,"normalpath(),sflag=1\n"); fflush(logp);} dsign=(VDENOM>=0.)?1:-1; /*set non-invert'y flags*/ if (v[0]*(ZDOT-v[0]/4.)-nu[0]*(z-zeta[0])*(z-zeta[0])>0.) eflag=12; else if (v[1]*(ZDOT-v[1]/4.)-nu[1]*(z-zeta[1])*(z-zeta[1])>0.) eflag=13; else { if (bdf!=2) d02cjc(neq, vprime, &z, v, s.zend, TOL, Nag_Mixed, out, gbothint, &comm, &fail); else d02ejc(neq, vprime, jacobboth, &z, v, s.zend, TOL, Nag_Mixed, out, gbothint, &comm, &fail); } } else if ((sflag==2)||(sflag==3)) /*1 player interior*/ { if(verbose==1) {fprintf(logp,"normalpath(),sflag=%d\n",sflag); fflush(logp);} dsign=(xi[sflag-2]+v[sflag-2]/2.>=beta*z)?1:-1; /*for non-invert'y check*/ whocorn=(sflag==2)?1:0; /*identifies cornered player for g1int()*/ if (bdf!=2) { if(verbose==1) {fprintf(logp,"calling d02cjc\n"); fflush(logp);} d02cjc(neq, vprime, &z, v, s.zend, TOL, Nag_Mixed, out, g1int, &comm, &fail); } else { if(verbose==1) {fprintf(logp,"calling d02ejc\n"); fflush(logp);} d02ejc(neq, vprime, jacob1, &z, v, s.zend, TOL, Nag_Mixed, out, g1int, &comm, &fail); } } else if (sflag==4) /*both cornered, no singularity possible*/ { if(verbose==1) {fprintf(logp,"normalpath(),sflag=4\n"); fflush(logp);} k[0]=(v[0]-2*nu[0]*(zeta[0]/(beta+delta[0])-z/delta[0]))/pow(z,(beta+delta[0])/beta); k[1]=(v[1]-2*nu[1]*(zeta[1]/(beta+delta[1])-z/delta[1]))/pow(z,(beta+delta[1])/beta); if((v[0]=2*nu[0]*(zeta[0]-z)/(beta+delta[0])) eflag=6; /*sufficient to stay cornered as per 15/10/99 version of jobmkt99.tex*/ /*n.b. corner checks may only work when increasing in z*/ else { ki=k[0]>=k[1]?k[1]:k[0]; whichi=(k[0]>=k[1])?1:0; c05adc(z+TOL,zend,&zevent,kroot,TOL,0.0,&fail); /*finds where 1 leaves corner*/ z=zevent; if (z>zend) eflag=6; else if (k[0]=k[1]) /*both leave together*/ { v[0]=v[1]=k[0]*pow(z,-(beta+delta[0])/beta)+2.*nu[0]*(zeta[0]/(beta+delta[0]) -z/(2.*beta+delta[0])); } else /*1 leaves first*/ { v[0]=k[0]*pow(z,-(beta+delta[0])/beta)+2.*nu[0]*(zeta[0]/(beta+delta[0]) -z/(2.*beta+delta[0])); v[1]=k[1]*pow(z,-(beta+delta[1])/beta)+2.*nu[1]*(zeta[1]/(beta+delta[1]) -z/(2.*beta+delta[1])); } } } if(verbose==1) {fprintf(logp,"eflag=%d,%s\n",eflag,fail.message); fflush(logp);} if (path==1) {fprintf(logp,"%8.2f %20.16f %20.16f\n",z,v[0],v[1]);fflush(logp);} } /*end of while loop*/ /*outputting results to result logfile*/ /* Structure: NE_NO_SIGN_CHANGE checked first as it would be caught by NE_NOERROR. * NE_NOERROR checked next as eflag might have been tripped before the error occurred. */ if ((fail.code==NE_NO_SIGN_CHANGE)&&(eflag==0)) (pcflag==0?fprintf(logr,"T,"):fprintf(logr,"t,")); else if (fail.code!=NE_NOERROR) { if (bdf==1) bdf=2; else if (bdf==2) {bdf=1; pcflag==0?(fprintf(logr,"F,")):(fprintf(logr,"f,"));} else {pcflag==0?(fprintf(logr,"F,")):(fprintf(logr,"f,"));} if(verbose==1) {fprintf(logp,"bdf=%d\n",bdf); fflush(logp);} } else { if (bdf==2) bdf=1; /*resets bdf in case d02ejc trips eflag but not NE_ TOL_PROGRESS*/ if (eflag==4) (pcflag==0?fprintf(logr,"S,"):fprintf(logr,"s,")); else if (eflag==6) (pcflag==0?(fprintf(logr,"C,")):(fprintf(logr,"c,"))); else if (eflag==10) (pcflag==0?(fprintf(logr,"E,")):(fprintf(logr,"e,"))); else if (eflag==11) (pcflag==0?(fprintf(logr,"Q,")):(fprintf(logr,"q,"))); else if (eflag==12) (pcflag==0?(fprintf(logr,"V,")):(fprintf(logr,"v,"))); else if (eflag==13) (pcflag==0?(fprintf(logr,"W,")):(fprintf(logr,"w,"))); else if (eflag==15) (pcflag==0?(fprintf(logr,"U,")):(fprintf(logr,"u,"))); else if (eflag==20) (pcflag==0?(fprintf(logr,"P,")):(fprintf(logr,"p,"))); fflush(logr); } if (bdf!=2) /*don't reset log files if haven't printed yet*/ { fflush(logr); if (path==1) {fprintf(logp,"%8.2f %20.16f %20.16f\n\n",z,v[0],v[1]); fflush(logp);} } /* Updating the state pointer before returning control to the calling routine * so that paths can be continued if desired. Maybe unnecessary in case of v[] * as names of arrays are automatically pointers? */ *state_ptr=z; *(state_ptr+1)=v[0]; *(state_ptr+2)=v[1]; } static double setpara() /********************************************************************************\ * This routine reads the parameter values, including the programme options, from * * a data file. * \********************************************************************************/ { FILE *in; static char input[40]; double zstart; if(verbose==1) {fprintf(logp,"setpara() entered\n"); fflush(logp);} if ((in=fopen("vf.data","r"))==NULL) { printf("Unable to open datafile\n."); exit(197); } else { fgets(input,40,in); beta=atof(input); fgets(input,40,in); delta[0]=atof(input); fgets(input,40,in); nu[0]=atof(input); fgets(input,40,in); xi[0]=atof(input); fgets(input,40,in); zeta[0]=atof(input); fgets(input,40,in); delta[1]=atof(input); fgets(input,40,in); nu[1]=atof(input); fgets(input,40,in); xi[1]=atof(input); fgets(input,40,in); zeta[1]=atof(input); fgets(input,40,in); zend=atof(input); fgets(input,40,in); uflag=atoi(input); fgets(input,40,in); initcond=atoi(input); fgets(input,40,in); int_step=atoi(input); fgets(input,40,in); ul=atof(input); fgets(input,40,in); ll=atof(input); fgets(input,40,in); ul2=atof(input); fgets(input,40,in); ll2=atof(input); fgets(input,40,in); zstart=atof(input); fgets(input,40,in); strcpy(scen_name,input); /*scenario name*/ if((delta[0]==delta[1])&&(nu[0]==nu[1])&&(xi[0]==xi[1])&&(zeta[0]==zeta[1])) symmetric=1; /*reads symmetry switch directly from the datafile*/ } fclose(in); return zstart; } static void out(Integer neq, double *zsol, double v[], Nag_User *comm) /********************************************************************************\ * Prints intermediate values from d02cjc() to logfile. Frequency as specified * * datafile. * \********************************************************************************/ { Integer i; double x[2]; struct user *s = (struct user *)comm->p; if(verbose==1) {fprintf(logp,"out() entered\n"); fflush(logp);} for (i=0;i=-2.*xi[i])?(xi[i]+v[i]/2.):0); /*eases utility plotting*/ if (path==1) { fprintf(logp,"%8.2f", *zsol); /*outputs z, v0, v1*/ for (i=0; izend - (double)s->k * s->h; s->k--; } static double gbothint(Integer neq, double z, double v[], Nag_User *comm) /********************************************************************************\ * Checks the path when both players play x_i > 0 against a series of events. * * Path conditioning is assessed independently of the eflag checks. * \********************************************************************************/ { int wflag; if(verbose==1) { fprintf(logp,"%8.2f %20.16f %20.16f, gbothint() entered\n",z,v[0],v[1]); fflush(logp); } if ((v[0]<-2.*xi[0])&&(v[1]<-2.*xi[1])) { if (((deriv(0,z,v,1))*(deriv(0,z,v,4))<0.)||((deriv(1,z,v,1))*(deriv(1,z,v,4))<0.)) wflag=gevent(11); /*quasi-non-invertibility check*/ else wflag=sflag=4; /*passes through corner w/o quasi-n-i*/ if(verbose==1) { fprintf(logp,"int: v'=(%8.6f,%8.6f), corn: v'=(%8.6f,%8.6f)\n", deriv(0,z,v,1),deriv(1,z,v,1),deriv(0,z,v,4),deriv(1,z,v,4)); fflush(logp); } } else if (v[0]<-2.*xi[0]) { /*only check quasi-n-i for cornering player: other can change slope w/o doubling back*/ if ((deriv(0,z,v,1))*(deriv(0,z,v,3))<0.) wflag=gevent(11); else wflag=sflag=3; if(verbose==1) { fprintf(logp,"int: v'=(%8.6f,%8.6f), corn: v'=(%8.6f,%8.6f)\n", deriv(0,z,v,1),deriv(1,z,v,1),deriv(0,z,v,3),deriv(1,z,v,3)); fflush(logp); } } else if (v[1]<-2.*xi[1]) { if ((deriv(1,z,v,1))*(deriv(1,z,v,2))<0.) wflag=gevent(11); /*quasi-n-i check*/ else wflag=sflag=2; if(verbose==1) { fprintf(logp,"int: v'=(%8.6f,%8.6f), corn: v'=(%8.6f,%8.6f)\n", deriv(0,z,v,1),deriv(1,z,v,1),deriv(0,z,v,2),deriv(1,z,v,2)); fflush(logp); } } else if ((v[0]*(ZDOT-v[0]/4.)-nu[0]*(z-zeta[0])*(z-zeta[0])>0.)&& (v[1]*(ZDOT-v[1]/4.)-nu[1]*(z-zeta[1])*(z-zeta[1])>0.)) wflag=gevent(15); else if (v[0]*(ZDOT-v[0]/4.)-nu[0]*(z-zeta[0])*(z-zeta[0])>0.) wflag=gevent(12); else if (v[1]*(ZDOT-v[1]/4.)-nu[1]*(z-zeta[1])*(z-zeta[1])>0.) wflag=gevent(13); else if (VDENOM*dsign<=0.) wflag=gevent(4); else if ((fabs(VDENOM)<=TOL)&&!reverse) wflag=gevent(10); /*small denominator in ode*/ else wflag=0; if (((VDENOM>((v[0]-v[1])*v[0]/4.))?((VDENOM>((v[1]-v[0])*v[1]/4.))? (VDENOM+v[0]*v[1]/4.):(v[1]*v[1]/4.)): ((fabs(v[0])>fabs(v[1]))?(v[0]*v[0]/4.):(v[1]*v[1]/4.)))>fabs(VDENOM)*TCOND) pcflag=1; /*poorly conditioned*/ if (((eflag==4)||(eflag==10))&&(fabs(ZDOT)<=TOL)) wflag=gevent(20); /*possible cusp*/ return (wflag==0)?1.:-1.; } static double g1int(Integer neq, double z, double v[], Nag_User *comm) /********************************************************************************\ * Checks the path when player i is cornered against a series of events. As sflag* * is altered by tentative steps of g1int() it cannot be used to track the phase * * in the interim. As in gbothint() the pcflag is now distinct from the eflag. * \********************************************************************************/ { Integer i,wflag; if(verbose==1) { fprintf(logp,"%8.2f %20.16f %20.16f, g1int() entered\n",z,v[0],v[1]); fprintf(logp,"z^{.}=%20.12f\n",ZDOT); fflush(logp); } i=whocorn; /*i is the cornered player*/ if (v[i]>=-2.*xi[i]) { /*only check quasi-n-i for cornering player: other can change slope w/o doubling back*/ if ((deriv(i,z,v,(i==1)?2:3))*(deriv(i,z,v,1))<0.) wflag=gevent(11); /*quasi-n-i check*/ else wflag=sflag=1; if(verbose==1) { fprintf(logp,"int: v'=(%8.6f,%8.6f), corn: v'=(%8.6f,%8.6f)\n", deriv(0,z,v,1),deriv(1,z,v,1),deriv(0,z,v,(i==1)?2:3),deriv(1,z,v,(i==1)?2:3)); fflush(logp); } } /*asymptotic approach to corner singularity checks*/ else if (dsign*ZDOT<=0.) wflag=gevent(4); /*denominator changes sign*/ else if ((fabs(ZDOT)<=TOL)&&(!reverse)) wflag=gevent(10); /*small ODE denominator*/ else if((-xi[i]*xi[i]-nu[i]*(z-zeta[i])*(z-zeta[i])+v[i]*ZDOT>0.)&& (v[1-i]*(ZDOT-v[1-i]/4.)-nu[1-i]*(z-zeta[1-i])*(z-zeta[1-i])>0.)) wflag=gevent(15); else if (-xi[i]*xi[i]-nu[i]*(z-zeta[i])*(z-zeta[i])+v[i]*ZDOT>0.) wflag=gevent((i==0)?12:13); else if (v[1-i]*(ZDOT-v[1-i]/4.)-nu[1-i]*(z-zeta[1-i])*(z-zeta[1-i])>0.) wflag=gevent((i==0)?13:12); else if (v[1-i]<-2.*xi[1-i]) wflag=4; /*IS THIS RIGHT?*/ else wflag=0; if ((v[i]*v[i]/(2.*ZDOT*ZDOT))>TCOND) pcflag=1; /*poor conditioning*/ return (wflag==0)?1.:-1.; } static int gevent(Integer code) /********************************************************************************\ * When an event has been triggered in gbothint or g1int this updates eflag. * \********************************************************************************/ { if(verbose==1) {fprintf(logp,"gevent() entered with code=%d\n",code); fflush(logp);} eflag=code; return code; /*UGLY: UNCLEAR WHAT wflag DOES; TOO MUCH INFO?*/ } static void vprime(Integer neq, double z, double v[], double vp[], Nag_User *comm) /********************************************************************************\ * The v' values, generally. Called only by d02cjc() * \********************************************************************************/ { int i; for (i=0;i-2.*xi[0])&&(v[1]>-2.*xi[1])) sflag=1; /*both interior*/ else if (v[0]>-2.*xi[0]) sflag=2; /*1 in interior*/ else if (v[1]>-2.*xi[1]) sflag=3; /*2 in interior*/ else sflag=4; } static int lpfiler(i,j,lastdot) /********************************************************************************\ * Generates the dots in .lp file from which results diagrams can be generated. * * The lastdot integer slims the .lp file by only printing when the dot differs * * from its predecessor. * \********************************************************************************/ { int thisdot; if ((eflag==4)||(eflag==10)) { if(lastdot!=2) {fprintf(loglp,"\\drawthickdot{%d}{%d}\n",i,j); fflush(loglp);} thisdot=2; } else if ((eflag==6)||(eflag==12)||(eflag==13)||(eflag==15)) { if(lastdot!=1) {fprintf(loglp,"\\drawdot{%d}{%d}\n",i,j); fflush(loglp);} thisdot=1; } else if (eflag==11) { if(lastdot!=0) {fprintf(loglp,"\\drawhollowdot{%d}{%d}\n",i,j); fflush(loglp);} thisdot=0; } else if (eflag==14) { if(lastdot!=3) {fprintf(loglp,"\\drawhollowdot{%d}{%d}\n",i,j); fflush(loglp);} thisdot=3; } else thisdot=-1; return thisdot; } static void jacob1(Integer neq, double z, double v[], double pw[], Nag_User *comm) /********************************************************************************\ * Provides Jacobian information on the system of ODEs when 1 player is cornered * * for use in the BDF routine d02ejc (e.g. for sflag == 2 || 3). * \********************************************************************************/ { #define PW(I,J) pw[((I)-1)*neq+(J)-1] double vjpvj, vjpvi, vipvj, vipvi; /*partial derivatives*/ int i=((sflag==2)?1:0); /*interior player*/ vjpvi = 0.; vjpvj = (beta+delta[1-i]-((beta+delta[1-i])*v[1-i]+2.*nu[1-i]*(z-zeta[1-i]))/(2.*ZDOT))/ZDOT; vipvi = (beta+delta[i]-((beta+delta[1-i])*v[1-i]+2.*nu[1-i]*(z-zeta[1-i]))/(2.*ZDOT))/ZDOT; vipvj = -((beta+delta[i])*v[i]+2.*nu[i]*(z-zeta[i])+beta+delta[1-i] -((beta+delta[1-i]*v[1-i]+2.*nu[1-i]*(z-zeta[1-i])))/ZDOT)/(2.*ZDOT*ZDOT); if (sflag==2) {PW(1,1)=vjpvj; PW(1,2)=vjpvi; PW(2,1)=vipvj;PW(2,2)=vipvi;} else if (sflag==3) {PW(1,1)=vipvi; PW(1,2)=vipvj; PW(2,1)=vjpvi;PW(2,2)=vjpvj;} } static void jacobboth(Integer neq, double z, double v[], double pw[], Nag_User *comm) /********************************************************************************\ * Provides Jacobian information on the system of ODEs when both players are * * cornered for use in the BDF routine d02ejc. * \********************************************************************************/ { #define PW(I,J) pw[((I)-1)*neq+(J)-1] int i; for (i=0;i