Shadowing Chaotic Motion

Stan Wagon

Macalester College, St. Paul, Minnesota

Rob Knapp

Wolfram Research, Inc., Champaign, Illinois

It pays to be imaginative when combining Mathematica's capabilities. Here optimization and high-precision are combined in a way that yields information about shadowing and the qualitative reliability of machine precision calculations.

Shadowing a Noisy Trajectory: The quadratic map

Consider the quadratic map r x (1-x) where r=4. 
Start at 0.1 and iterate using machine precision: error quickly builds up as the roundoff propagates.
Here blue is a true trajectory from 1/10 using high-precision. 
Red is the noisy trajectury; it is not a true trajectory. Or is it?

f[x_] :=4x (1-x);
ListLinePlot[{NestList[f,N[1/10,100],100],NestList[f,0.1,100]}, Frame->True, Axes->False, AspectRatio-> 0.4,FrameTicks->{Automatic, Automatic,None,None}, PlotStyle->{{Blue,Thickness[0.003]},{Lighter@Red,Thickness[0.007]}}, ImageSize->500, BaseStyle->14,PlotLabel->"Red = machine precision from 0.1, blue = high precision from 1/10"]

The shadowing question: Is there is an initial value whose true trajectory is essentially identical to the noisy trajectory. We will describe a new way of answering such questions. The point is that when this happens, then this means the noisy trajectory is, within a small tolerance, the same as a true trajectory, but of another starting value.

This demo gets the idea across. One can, slowly but surely, use high precision to zero in on the shadow seed.

Module[{NN, makeTimes, f, trajectory, base, baseplot,r},
NN[xx_] := SetPrecision[xx,80];
limit4=-1.11566519693586352887390764612 10^-17;
steps = 
{-.111, -.11156, -.1115665,-.111566519,-.11156651969,-.1115665196935,-.111566519693586,-.11156651969358635};
makeTimes[z_]:=Style[z,FontFamily->"Times"];
f[r_][xx_]:=r xx(1-xx);
trajectory[r_,seed_]:=NestList[f[r],seed,100];
base[s_]:=base[s]=trajectory[s,0.1];
baseplot[z_]:=ListLinePlot[base[z],PlotStyle->{Thickness[0.016],Red}];
olddiff=0;
Manipulate[
If[update,olddiff+=NN[change (10^scale)];change=0;update=False];
If[autosteps>0, change=steps[[autosteps]];scale=-16];
new=trajectory[r,NN[1/10+olddiff]+NN[change (10^scale)]];
Pane[Column[{Row[{
makeTimes[StringForm["Difference from (1/10):  ``\n",Style[ToString[SetPrecision[olddiff+change (10^scale),30],TraditionalForm]]]],
makeTimes[StringForm["Total error as (Sqrt["sum of squares"]):  ``\n",Style[ToString[NumberForm[Norm[base[r]-new],8],TraditionalForm]]]]}],
Show[baseplot[r],
ListLinePlot[new,PlotStyle->{Thickness[0.005],Blue}],
Frame->True,Axes->False,AspectRatio->0.4,ImageSize->500,FrameTicks->{Automatic,{0,0.5,1},None,None},BaseStyle->14,PlotLabel->None,PlotRange->{{start,Min[start+len,100]},{0,1}}]}],ImageSize->{510,310}],{{scale,-16,"power of 10 scale for change"},Range[-8,-30,-2]},{{change,0.,"change (without power of 10)"},-1,1,Appearance->"Closed"},{{start,40,"display start"},1,99,1},{{len,40,"display length"},1,100,1},{update,{False,True}},{{r,4,"r"},{4,38/10->3.8}},{{autosteps,0,"Automatic convergence steps"}, {0,1,2,3,4,5,6,7,8}, SetterBar},TrackedSymbols:>{scale,r,update,start,len,change, autosteps}, SaveDefinitions->True]]

There is nothing surprising here because of a certain bit-shift interpretation that is known to hold for 4 x (1-x). But everything holds for other cases as well, such as (38/10) x (1-x).

Shadowing a Noisy Trajectory: The quadratic map

Here is a reminder that the chaotic values of r are interspersed with more stable values.

Clear[r];
Options[BifurcationPlot]={PlotPoints->2000,PlotStyle->{Black,PointSize[0.001]}};
BifurcationPlot[f_,{r_,a_,b_},{x_,x0_},{iter0_,iterShow_},opts___]:=Module[{sty,n,makePts,rVals,data},
{sty,n}={PlotStyle,PlotPoints}/.{opts}/.Options[BifurcationPlot];
makePts[{s_,v_}]:=({s,#}&)/@v;
cf=Compile[{{r,_Real,1},{x,_Real,1}},Evaluate[f]];
rVals=Range[N[a],b,((b-a)/(n-1))];
data=Transpose[{rVals,Transpose[NestList[cf[rVals,#]&,Nest[cf[rVals,#1]&,Array[x0&,n],iter0], iterShow]]}];
Graphics[Append[Flatten[{sty}],Point[Flatten[makePts/@data,1]]],Sequence@@FilterRules[{opts},Options[Graphics]], 
AspectRatio->1/3,Frame->True,FrameTicks->{Automatic,Range[0,1,.5], None, None},Axes->None]];
BifurcationPlot[r x(1-x),{r,2.95,4},{x,0.5},{500,100},FrameTicks->{{3,3.45,3.57,3.8,4},Automatic,{},None},GridLines->{{3,3.45,3.57,3.8,3.839},None}, BaseStyle->16]

The point: It is in many cases known that shadowing exists. In other cases it is not known. We describe here techniques using optimization that can explicitly find the shadows.

2. Automating the Optimization

r=4; fQuad[x_]:=r  x (1-x);
der[x_]:=Evaluate[fQuad'[x]];
fIter[x_,n_]:=Nest[fQuad,x,n];
fTraj[x_,n_]:=NestList[fQuad,x,n];

The loss of precision is 60 digits for 100 iterations (this is actually too pessimistic as the error is overestimated).

Precision[fIter[N[(1/10),500],100]]

fIter[0.1,100]
fIter[N[(1/10),65],100]
val100=fIter[N[(1/10),100],100]

A key step in optimization is understanding the derivative of the objective. That is very easy to compute here.

sensitivity[x_,n_]:=Most[FoldList[Times,{1},der[fTraj[x,n]]]]

ListPlot[Re[Log[10,Flatten[sensitivity[N[(1/10),70],100]]]],Frame->True,BaseStyle->14,FrameTicks->{Automatic,Range[0,25,5],{},{}},PlotStyle->{PointSize[0.015],Thickness[0.01]},Mesh->All,MeshStyle->{PointSize[0.008],Yellow},FrameLabel->{"Iterations","(Subscript[log, 10]) sensitivity"}]

Reality check:

sensitivity[N[(1/10),70],100][[-1]]
ϵ=(10^(-40));(1/ϵ)(fIter[N[(1/10)+ϵ,100],100]-val100)

2. Automating the Optimization

Get noisy trajectory.

noisy=SetPrecision[fTraj[0.1,100],100];  noisy[[-1]]

Set up the search:  

• Define a residual function to give the vector of residuals when comparing noisy to true. Requires upgrading the precision of the noisy trajectory to match the working precision. WP can be taken to be number of iterations.

• Use FindMinimum where the objective is specified in the Method option: Levenberg-Marquardt is made for minimizing norms of residual vectors (so the 0 as argument is irrelevant (weird))

• We also specify the sensitivity at each iteration by setting of the Jacobian option.

Shadow[noisy_, n_, start_] := (
residual[x0_?NumericQ] :=Block[{$MinPrecision = n,$MaxPrecision = n},
fTraj[x0,n]-Take[noisy,n+1]];
FindMinimum[0,{x0,SetPrecision[start, n]}, WorkingPrecision->n,
Method->{"LevenbergMarquardt","Residual"->residual[x0],"Jacobian"->sensitivity[x0,n]}])

Preliminary: How Much to Bite Off?

Key point: Optimization works if the objective is well-behaved (quadratic) around the answer. If one tries to shadow for 100 iterations, the residual is not at all well-behaved. Here is the neighborhood of the answer at a scale of (10^(-27)).

residual1[x0_?NumericQ,n_]:=fTraj[x0,n]-Take[noisy,n+1];
sol50=Shadow[noisy,50,1/10];
shadowSeed50=x0/.sol50[[2]];
Do[sol=Shadow[noisy,nn,If[nn==50,1/10,x0/.sol[[2]]]],{nn,50,100,50}];
shadowSeed100=x0/.sol[[2]];
sol


length=100;p=27;npts=32000;
data100=Table[Max[Abs[residual1[N[x0,length],length]]],{x0,shadowSeed100-(10^(-p)),shadowSeed100+(10^(-p)),2 (10^(-p))/npts}];
ListPlot[data100,Frame->True,PlotStyle->PointSize[0.003],Axes->False,PlotRange->{-0.1,1.1},FrameTicks->{{{npts/2,"shadow seed"},{0,StringForm["seed - (10^(-``))",p]},{npts,StringForm["seed + (10^(-``))",p]}},{0,0.5,1},None,None},Epilog->{PointSize[0.02],Red,Point[{npts/2,Max[Abs[residual1[shadowSeed100,100]]]}]}, BaseStyle->14]

Since we have no idea where the seed is, this is impossible to find. The shadow seed is about (10^(-17)) left of (1/10). Here it is at a scale of (10^(-16)) where we can see both (1/10) and the shadow seed. The answer is a needle in a giant haystack. Starting at (1/10), one has no way to find the shadow seed.

length=100;p=16;npts=32000/4;
data100Large=Table[Max[Abs[residual1[N[x0,length],length]]],{x0,shadowSeed100-(10^(-p)),shadowSeed100+(10^(-p)),2 (10^(-p))/npts}];
ListPlot[data100Large,Frame->True,PlotStyle->PointSize[0.003],Axes->False,PlotRange->{{0, npts .6},{-0.1,1.1}},FrameTicks->{{{npts/2,"shadow\nseed"},{0,StringForm["seed - (10^(-``))",p]},{npts,StringForm["seed + (10^(-``))",p]},
{ (npts/2)+((1.11/10) npts/2), 1/10}},{0,0.5,1},{{ (npts/2)+((1.11/10) npts/2), 1/10}},None},Epilog->{PointSize[0.02],Red,Point[{npts/2,Max[Abs[residual1[shadowSeed100,100]]]}]}, BaseStyle->14, GridLines->{{(npts/2)+((1.11/10) npts/2)},{}}]

Aside: How Much to Bite Off?

But if we cut back to only 50 iterations then the residual is not so wild, and it is easy for FindMinimum to do its job. Here it is at (10^(-15)) scale.

p=15;npts=500;length=50;
data=Table[Max[Abs[residual1[N[x0,length],length]]],{x0,shadowSeed50-(10^(-p)),shadowSeed50+(10^(-p)),2 (10^(-p))/npts}];
ListPlot[data,Frame->True,PlotStyle->PointSize[0.003],Axes->False,PlotRange->{0, 1},FrameTicks->{{{npts/2,"shadow seed"},{0,StringForm["seed - (10^(-``))",p]},{npts,StringForm["seed + (10^(-``))",p]}},{0,0.5,1},None,None},Epilog->{PointSize[0.03],Red,Point[{npts/2,Max[Abs[residual1[shadowSeed50,length]]]}]}, PlotRangePadding->{10,.05}, BaseStyle->14]

And here is a scale that shows the answer and the start of the search, (1/10)

p=16;npts=500;length=50;
data=Table[Max[Abs[residual1[N[x0,length],length]]],{x0,shadowSeed50-(10^(-p)),shadowSeed50+(10^(-p)),2 (10^(-p))/npts}];
ListPlot[data,Frame->True,PlotStyle->PointSize[0.003],Axes->False,PlotRange->{{0, npts .6},{0, .3}},FrameTicks->{{{npts/2,"shadow\nseed"},{0,StringForm["seed - (10^(-``))",p]},{npts,StringForm["seed + (10^(-``))",p]}(*,{ (npts/2)+((1.11/10) npts/2), 1/10}*)},{0,0.3},{{ (npts/2)+((1.11/10) npts/2), 1/10}},None},Epilog->{PointSize[0.03],Red,Point[{npts/2,Max[Abs[residual1[shadowSeed50,length]]]}]}, PlotRangePadding->{10,.05}, BaseStyle->10, GridLines->{{(npts/2)+((1.11/10) npts/2)},{}}]

So the idea is to bootstrap up from 50 in steps of 50.

2. Automating the Optimization

sol50=Shadow[noisy,50,1/10]

shadowSeed50=x0/.sol50[[2]]

And now we iterate as the number of steps starts at 50, and then jumps to 100.

Do[sol=Shadow[noisy,nn,If[nn==50,1/10,x0/.sol[[2]]]],{nn,50,100,50}]
shadowSeed100=x0/.sol[[2]]

shadowTraj=SetPrecision[fTraj[shadowSeed100,100],100];
Max[Abs[shadowTraj-noisy]]

Here the noisy trajectory is red; the high precision is the black line. In fact, it is visually good out to 150.

Show[ListLinePlot[fTraj[0.1,165], PlotStyle->{Thickness[0.01],Red}],
ListLinePlot[fTraj[shadowSeed100,165], PlotStyle->{Thickness[0.004], Green}], Frame->True, AspectRatio->0.4, PlotRange-> {-0.1, 1.1}, Axes->False, FrameTicks->{Automatic, {0,0.5,1}, None, None}, GridLines->{{150},{}}, BaseStyle->14, ImageSize->500]

There is no trouble using this technique to get a shadow value for 500 or 1000 iterations. One can ask questions such as: Is the shadow seed equally likely to be left as right?

All this works the same way when r is, say, 38/10.

3. The Duffing equation

A classic sensitive differential equation is the Duffing equation, a forced oscillator.

Clear[t,x]; duff={x''[t]+(15/100) x'[t]-x[t]+((x[t])^3)==(3/10) Cos[t]};

T=100;xStart=0;xpStart=0;  
basicSol=x[t]/.First[NDSolve[{duff,x[0]==xStart,x'[0]==xpStart},x,{t,T}, InterpolationOrder->All]]

ParametricPlot[Evaluate[{basicSol,D[basicSol,t]}], {t,0,T}, Frame->True, PlotStyle->{Thick, Red}, Epilog->{PointSize[Large], Point[{xStart, xpStart}]}, FrameTicks->{Range[-1,1,1],Range[-1,1,1],None,None},BaseStyle->14]

3. The Duffing equation

Here is the high-precision solution from the exact initial value (0,0).

highPrecSol=x[t]/.First[
NDSolve[{duff,x[0]==xStart,x'[0]==xpStart},x,{t,T}, InterpolationOrder->All,PrecisionGoal->10, WorkingPrecision->100]];

Comparing:

ParametricPlot[Evaluate[{{basicSol,D[basicSol,t]},{highPrecSol,D[highPrecSol,t]}}],{t,0,T},Frame->True,PlotStyle->{{Thickness[0.02],Red},{Thickness[0.005],Blue}},Axes->False, BaseStyle->14,Epilog->{PointSize[0.03],Point[{xStart,xpStart}]},FrameTicks->{{-1,0,1},{-1,0,1},None,None}]

Or look just at x and not the derivative.

Plot[{basicSol,highPrecSol},{t,0,T},Frame->True,PlotStyle->{{Thickness[0.02],Red},{Thickness[0.005],Blue}},Axes->False,FrameTicks->{Range[0,100,20],{-1,0,1},None,None}, BaseStyle->14]

3. The Duffing equation

Now, a key step is working out the sensitivity of the solution to changes in the initial value. This can then be used as a Jacobian when we use LevenbergMarquardt to optimize. Rob worked this out: the info we want is inherent in the DE, which one considers as a function of t and also the parameters x0 and xp0.

General method for sensitivity

Write the equations in terms of not only the usual t, but all of the variables you want sensitivity for:

ivar=t;
origdvars={x};
origeqns={x''[t]+((15 x'[t])/100)-x[t]+((x[t])^3)==((3 Cos[t])/10),x[0]==x0,x'[0]==xp0};
params={x0,xp0};

Arguments of the modified functions will be f[p1, p2, ..., t], so the sensitivity is given by derivatives with respect to the parameter arguments.

parg[iv_] := Apply[Sequence, Append[params, iv]];
args = parg[ivar];
dvars = Map[#[args]&,{x}];
sdvars = Flatten[{dvars, D[dvars, {params}]}]

Replace f[t] with f[p1,p2,…,t] and correspondingly for derivatives:

rules = Flatten[Map[
{#[iv_]:> #[parg[iv]], Derivative[j_][#][iv_]:>Apply[Derivative, Append[ConstantArray[0, Length[params]],j]][#][parg[iv]]}&, origdvars]]

eqns = origeqns /. rules

The senstivity equations are just the derivatives of the substituted equation with respect to the parameters

seqns=D[eqns,{params}]

Now just solve with NDSolve (for particular values of the parameters!!!)

Block[{x0=1,xp0=0},sol=NDSolve[{eqns,seqns},sdvars,{t,0,100}]]

Duffing equation with sensitivity equations for

       ((∂x(t))/(∂(Subscript[x, 0]))) represented by dx0                   ((∂x(t))/(∂(Subscript[x, 0]^′))) represented by dxp0

duff={x''[t]+(15/100) x'[t]-x[t]+((x[t])^3)==(3/10) Cos[t]};
sx0={dx0''[t]+(15/100) dx0'[t]-dx0[t] +3 ((x[t])^2)dx0[t] ==0,dx0[0]==1,dx0'[0]==0};
sxp0={dxp0''[t]+(15/100) dxp0'[t]-dxp0[t] (1-3 ((x[t])^2))==0,dxp0[0]==0,dxp0'[0]==1};

Here dx0[t] represents the rate of change of x[t] when one changes the initial value x(0). We know this equation is sensitive so we expect this rate of change to get quite large when t is, say 100. It turns out that for such t the sensitivity is about (10^8): So a change of (10^(-8)) in the initial value can cause a change of 1 in x(100). The trick will be to solve these three equations simultaneously, thus getting good info about this sensitivity.

3. The Duffing equation

Time, sampling interval

T=100;Δt=1;samples=Range[0,T,Δt];

This means we are aiming for a match up to time 100 and we will sample just every unit to try to get the match. This appears to be good enough, though it would be no trouble to reduce Δt.

Working Precision, AccuracyGoal, PrecisionGoal for NDSolve

ag and pg should be less than half of wp to be sure roundoff error is not an issue. While 17 might be enough for this to work, the extra time in using 20 or even 50 is quite small. So let us be safe and use 30.

wp=30;ag=pg=10;xStart=0;xpStart=0;

The Basic Solution

Here is the basic noisy machine precision solution at the sample points only. The SetPrecision command is used to turn the low-precision numbers into wp-precision numbers for later use with FindMinimum. The InterpolationOrder->All setting changes how NDSolve saves its solution so that even between steps the solution is just as precise. This is needed here since the low and high precision solutions will use very different step sizes.

basicSol=x[t]/.First[NDSolve[{duff,x[0]==xStart,x'[0]==xpStart},x,{t,T}, InterpolationOrder->All]]

base = SetPrecision[Transpose[
First[{x[samples], x'[samples]} /. NDSolve[{duff,x[0]==xStart,x'[0]==xpStart},x,{t,T}, InterpolationOrder->All]]], wp];
base//Length

basePts=Transpose[{samples,First/@base}];
Show[Plot[basicSol,{t,0,T},PlotStyle->Thick],ListPlot[basePts,PlotStyle->{Red,PointSize[0.015]}],Frame->True,BaseStyle->14, Axes->False, FrameTicks->{Range[0,100,20],{-1,0,1}, {},{}}]

This shows that 100 sampling points comes close to capturing the curve. Could redo all this with 200 sample points.

High precision sensitivity solution

We now get a high precision sensitivity solution -- cached since it's used for the Jacobian too. This is a soln to Duffing together with the eqns that give the sensitivity.

The Extrapolation method is a good choice here because it is an arbitrary order method.  If you have a method of order p, typically the step size that NDSolve is able to take scales roughly like s=(((c (10^(-t))))^((1/p))) where t  is a measure of the tolerance (you can effectively think of it as the setting for PrecisionGoal). This means that the number of steps is roughly proportional to (10^(t/p)). If you have a fixed order method (e.g. the classical fourth order RungeKutta method with p =4) then for t=10, you might use 300 steps, but with tolerance t=20 you would use 100000 steps.  On the other hand, if you can increase p dynamically as the tolerance increases the number of steps can be reduced substantially.

Clear[hpssol];
hpssol[{x0_,xp0_}]:=hpssol[{x0,xp0}]=
NDSolve[{duff,x[0]==SetPrecision[x0,wp],x'[0]==SetPrecision[xp0,wp],sx0,sxp0},{x,dx0,dxp0},{t,T},Method->"Extrapolation",InterpolationOrder->All,WorkingPrecision->wp,AccuracyGoal->ag,PrecisionGoal->pg][[1]]

Again, note the interplay here: we are solving not just for the true value of x, but simultaneously for the sensitivity analysis. Here we compute one. These take a few seconds.

hpssol[{0,0}]

Here is the high-precision sensitivity for this particular start value. It is qualitatively similar to the low-precision version. Logs used so the dips are just the zeros.

Plot[Log[10,Abs[dx0[t]]]/.hpssol[{0,0}],{t,0,T},PlotStyle->Thick,Frame->True,BaseStyle->14,Axes->False,FrameTicks->{Range[0,100,20],({#1,(("10")^(#))}&)/@Range[0,8,2],{},{}}]

As a check we perturb a little — by (10^(-10)), the order of the PrecisionGoal — and find that the solution is the same. This indicates that our WP, PG, AG are good enough for truth, at least up to time 100.

Show[Plot[Evaluate[x[t]/.hpssol[{xStart,xpStart}]],{t,0,100},PlotStyle->{Thickness[0.015],Red}],Plot[Evaluate[x[t]/.hpssol[{xStart,xpStart}+(1/(10^10))]],{t,0,100},PlotStyle->{Thickness[0.006],Yellow}], BaseStyle->14]

Note: The sensitivity matrix changes as the initial value changes. But hpssol is recomputed each time!

The Plan

Now the plan is to try to find initial conditions for which the true solution matches VERY CLOSELY the base solution. This shows that the numerical error in the base soln is not just garbage, but that it represents a real solution to the very same equation, but with different initial values. This phenomenon is called shadowing. It is suspected to hold for Duffing, but this is not proved. It is proved for some other systems. So the plan is to define an objective function that, for starting values x0 and xp0, computes the true high-precision solution and compares it with the base solution by computing a residual based on 100 sample points. Then the idea is to minimize the residual. There are some sophisticated minimization algorithms, especially when the residual is a sum of squares (Levenberg-Marquardt). But some derivative data is needed -- the Jacobian.

Residual

The following function, for given x0 and xp0, uses high precision to determine the total error when compared to the base solution -- this error is called the residual. Our goal is to find values so that residual is very small.

Clear[r];
r[{x0_?NumberQ,xp0_?NumberQ}]:=Flatten[Transpose[{x[samples],x'[samples]}/.hpssol[{x0,xp0}]]-base]

At the origin there is a large residual.

Norm[r[{0,0}]]

Pretty large. How should one change the values? It seems conceivable that the sensitivity info could help us decide how to move ... but there are 100 sample points, and only two points to move.

The Key Step: The Jacobian Function

FindMinimum works, but we need to use either LevenbergMarquardt or Newton's Method, feeding the Jacobian info -- the sensitivity -- to the method. J gives the Jacobian info. LM is faster than Newton. Rob will address.

J[{x0_?NumberQ,xp0_?NumberQ}]:=Transpose[
{Flatten[Transpose[{dx0[samples],dx0'[samples]}/.hpssol[{x0,xp0}]]],Flatten[Transpose[{dxp0[samples],dxp0'[samples]}/.hpssol[{x0,xp0}]]]}]

J[{0,0}][[{1,2,3,4}]]

This encapsulates the sensitivity info at the sample points.

Length[J[{0,0}]]
J[{0,0}][[2]]
J[{0,0}][[202]]
Length[samples]

So J consists of 101 pairs.

The Objective function

obj[{x0_?NumberQ, xp0_?NumberQ}] := Module[{res = r[{x0,xp0}]},(res.res)/2]

obj[{xStart,xpStart}]

You can clearly see deviation in the plot -- the picture underlying the nonzero objective at (0,0).

Show[Plot[Evaluate[
x[t]/.First[NDSolve[{duff,x[0]==xStart,x'[0]==xpStart},x,{t,T},InterpolationOrder->All]]], {t,0,T}, PlotStyle->Thickness[.013]],
Plot[Evaluate[{x[t]} /. hpssol[{0,0}]],{ t,0,T}, PlotStyle->{Thickness[.005],Pink}], BaseStyle->14, Frame->True, Axes->False, FrameTicks->{Range[0,100,20], {-1,0,1}, {},{}}]

Summary: The red is truth. The blue has numerical error that makes it wrong past 80. Note that things are a little different for other start values. For (-1,1) the divergence starts at 60. Here is the phase plane view.

Show[ParametricPlot[Evaluate[{x[t],x'[t]}/.hpssol[{0,0}]],{t,0,T},PlotStyle->Thick],ListPlot[base,PlotStyle->{Red,PointSize[0.012]}],Frame->True,PlotRangePadding->0.1,Axes->False, BaseStyle->14,FrameTicks->{{-1,0,1},{-1,0,1},None,None}]

Success

So, we now put it all together: Takes 33 seconds.

Timing[{min,minx}=FindMinimum[0,{{x0,xStart},{xp0,xpStart}},Method->{"LevenbergMarquardt","Residual"->r[{x0,xp0}],"Jacobian"->J[{x0,xp0}]},WorkingPrecision->wp,AccuracyGoal->0.8 ag,PrecisionGoal->0.8 pg,EvaluationMonitor:>Print[NumberForm[{obj[{x0,xp0}],x0,xp0},3]]]]

The shadow is close to the original seed:    2 (10^(-8)).
The residual has nicely small norm of (10^(-11)).
And visually the solutions are wonderfully consistent.

Show[ParametricPlot[Evaluate[{x[t],x'[t]} /. hpssol[{x0, xp0} /.minx]],{ t,0,T},
PlotStyle->{{Blue,Thick}}],ListPlot[base,PlotStyle->{Red, PointSize[0.012]}], Frame->True,Axes->False, FrameTicks->{Range[-1,1,1],Range[-1,1,1],None,None}, PlotRangePadding->0.1, BaseStyle->14]

Show[Plot[Evaluate[{x[t]}/.hpssol[{x0,xp0}/.minx]],{t,0,T},PlotStyle->{Blue,Thickness[.02]}],
Plot[basicSol,{t,0,T},PlotStyle->{Yellow,Thickness[0.006]}],Frame->True,Axes->False,BaseStyle->14,FrameTicks->{Range[0,100,20],{-1,0,1},{},{}}]

So we have a solution that shadows the computed solution, with initial value (shadow seed) at distance (10^(-8)) from the noisy source and total residual about (10^(-12)).

Harder Initial Values (-1,1)

For more sensitive initial values we need to bootstrap. It works if we go to 75 first and then continue. Let's use 50 to be consistent with quadratic map.

Start at 50:

duff={x''[t]+(15/100) x'[t]-x[t]+((x[t])^3)==(3/10) Cos[t]};
sx0={dx0''[t]+(15/100) dx0'[t]-dx0[t] +3 ((x[t])^2)dx0[t] ==0,dx0[0]==1,dx0'[0]==0};
sxp0={dxp0''[t]+(15/100) dxp0'[t]-dxp0[t] (1-3 ((x[t])^2))==0,dxp0[0]==0,dxp0'[0]==1};

r[{x0_?NumberQ,xp0_?NumberQ}]:=Flatten[Transpose[{x[samples],x'[samples]}/.hpssol[{x0,xp0}]]-base];

obj[{x0_?NumberQ, xp0_?NumberQ}] := Module[{res = r[{x0,xp0}]},(res.res)/2];
J[{x0_?NumberQ,xp0_?NumberQ}]:=Transpose[
{Flatten[Transpose[{dx0[samples],dx0'[samples]}/.hpssol[{x0,xp0}]]],Flatten[Transpose[{dxp0[samples],dxp0'[samples]}/.hpssol[{x0,xp0}]]]}];

Clear[hpssol];
hpssol[{x0_,xp0_}]:=hpssol[{x0,xp0}]=
NDSolve[{duff,x[0]==SetPrecision[x0,wp],x'[0]==SetPrecision[xp0,wp],sx0,sxp0},{x,dx0,dxp0},{t,T},Method->"Extrapolation",InterpolationOrder->All,WorkingPrecision->wp,AccuracyGoal->ag,PrecisionGoal->pg][[1]];

Here we go out to 50.

wp = 30;ag = pg = 10;xStart=-1;xpStart=1;
T=50;Δt=1;samples=Range[0,T,Δt];
basicSol=x[t]/.First[NDSolve[{duff,x[0]==xStart,x'[0]==xpStart},x,{t,T}, InterpolationOrder->All]]
base = SetPrecision[Transpose[
First[{x[samples], x'[samples]} /. NDSolve[{duff,x[0]==xStart,x'[0]==xpStart},x,{t,T}, InterpolationOrder->All]]], wp];
Norm[r[{0,0}]]
J[{-1,1}][[{1,2,3,4}]]
Timing[{min, minx} = FindMinimum[0, {{x0, xStart},{xp0, xpStart}}, Method->{"LevenbergMarquardt","Residual"->r[{x0,xp0}],"Jacobian"->J[{x0, xp0}]}, WorkingPrecision->wp, AccuracyGoal->0.8 ag, PrecisionGoal->0.8 pg, EvaluationMonitor:>Print[NumberForm[{obj[{x0,xp0}],x0, xp0},3]]] ]

Show[
Plot[Evaluate[{x[t]} /. hpssol[{x0, xp0} /.minx]],{ t,0,T},PlotStyle->{Blue,Thickness[.02]}],
Plot[basicSol,{ t,0,T}, PlotStyle->{Yellow,Thickness[.006]}],Frame->True, Axes->False, BaseStyle->14, FrameTicks->{Range[0,100,20], {-1,0,1},{},{}}]

And now we use the shadow seeds learned for time 50 — {x0 -> -0.999999967671633359656003680726, 
 xp0 -> 0.999999950363271277469162306091} — as initial values for the bootstrap step to time 100.

Clear[hpssol];
hpssol[{x0_,xp0_}]:=hpssol[{x0,xp0}]=
NDSolve[{duff,x[0]==SetPrecision[x0,wp],x'[0]==SetPrecision[xp0,wp],sx0,sxp0},{x,dx0,dxp0},{t,T},Method->"Extrapolation",InterpolationOrder->All,WorkingPrecision->wp,AccuracyGoal->ag,PrecisionGoal->pg][[1]];
T=100;samples=Range[0,T,Δt];
basicSol=x[t]/.First[NDSolve[{duff,x[0]==xStart,x'[0]==xpStart},x,{t,T}, InterpolationOrder->All]];
base = SetPrecision[Transpose[
First[{x[samples], x'[samples]} /. NDSolve[{duff,x[0]==xStart,x'[0]==xpStart},x,{t,T}, InterpolationOrder->All]]], wp];

Timing[{min, minx} = FindMinimum[0, {{x0, x0/. minx},{xp0, xp0/. minx}}, Method->{"LevenbergMarquardt","Residual"->r[{x0,xp0}],"Jacobian"->J[{x0, xp0}]}, WorkingPrecision->wp, AccuracyGoal->0.8 ag, PrecisionGoal->0.8 pg, EvaluationMonitor:>Print[NumberForm[{obj[{x0,xp0}],x0, xp0},3]]] ]

Show[ParametricPlot[Evaluate[{x[t],x'[t]}/.hpssol[{x0,xp0}/.minx]],{t,0,T},PlotStyle->{Blue,Thickness[0.02]}],ParametricPlot[Evaluate[{basicSol,(∂_t)basicSol}],{t,0,T},PlotStyle->{Yellow,Thickness[.006]}],
Graphics[{PointSize[Large], Point[{-1,1}]}],Frame->True,Axes->False,BaseStyle->14,FrameTicks->{Range[0,100,20],{-1,0,1},{},{}}]

Show[
Plot[Evaluate[{x[t]} /. hpssol[{x0, xp0} /.minx]],{ t,0,T},PlotStyle->{Blue,Thickness[0.02]}],
Plot[basicSol,{ t,0,T}, PlotStyle->{Yellow,Thickness[.006]}],Frame->True, Axes->False, BaseStyle->14, FrameTicks->{Range[0,100,20], {-1,0,1},{},{}}]

3. The SIAM Billiard Problem

From the SIAM 100-Digit Challenge: 

A particle bounces of mirrors of radius (1/3) centered at integer lattice points, and starts at ((1/2),(1/10)) traveling east at speed 1, where is it at time 10?

Solving it: Just use high-precision.
Machine precision: ok for a bit, then roundoff error makes it wrong.

base=PhotonTrajectory[{0.5,0.1},{1.,0},200];
Length[base]
pr =PlotRange[ListLinePlot[base, PlotRange->All]]+{{-1,1}, {-1,1}};  (* buggy with just ListPlot *) 
mirrors={RGBColor[1,1,.6],EdgeForm[Thin],Table[Disk[{i,j},1/3], {i,-20,20}, {j,-20,20}]};
baseimage=Graphics[{mirrors,Thickness[0.004],Line[base], {Red, PointSize[Large], Point[base[[{1,2,5,14}]]]}},Frame->True,PlotRange->{{-2,1},{-2,2}},PlotRangeClipping->True,BaseStyle->14,ImageSize->300,FrameTicks->{Range[-3,3],Range[-3,3],{0.5},{0.1}}]

In fact, machine precision, starting from (0.5, 0.1) is enough to get 5 digits of the answer (but not 10).

trajHP=PhotonTrajectory[N[{(1/2),(1/10)},300],{1,0},41,WorkingPrecision->200];
Manipulate[Graphics[{{
Pink,Thickness[0.013],Line[Take[trajHP,n]]},
{mirrors,Thickness[0.004],Line[Take[base,n]], {Red, PointSize[Large], Point[base[[{1,2,5,14}]]]}}},PlotRange->{{-1.5,1.75},{-1.9,1.9}}, ImageSize->250], {n,1,41,1, ControlPlacement->Left, ControlType->VerticalSlider}, TrackedSymbols->{n}, SaveDefinitions->True]

Shadowing: Is the noisy trajectory close to a true trajectory? Yes.

I made a demo allowing manual search for the shadow. Makes it almost obvious.

In fact the minimization of the residual can be automated, a key step being the computation of the sensitivity. 
A student of mine (Henrik Hakonsen) has carried out the details and it works. The state space is 4-dimensional (position, direction). The tricky part is getting the Jacobian which requires multiplying matrices obtained from solving the quadratic equations of the intersections.

Key idea: From the target noisy trajectory we know the sequence of circles and so can use algebra and the chain rule to get the various derivatives. Here is the result, using 500 digits of working precision in the optimization. The error in the two trajectories is about (10^(-15)).

seed[150]

Here is a look at how the error behaves at a scale of (10^(-17)), looking at error at 30 bounces.

Code

Block[{$MinPrecision=wp, $MaxPrecision=wp},
Clear[residual];
residual[y0_?NumericQ]:=(mPrint["res",y0];Flatten[(#[[2]]& /@ (First/@fBillTraj[{{1/2,y0},{1,0}},n,Max[n,wp]]))-
(#[[2]]& /@ First/@SetPrecision[Take[noisy,n+1],wp])])];

n=30;
N[Norm[residual[seed[150]]]]

n=20;yShadow=seed[150];

pow=17;
data=Table[{-(y-(yShadow+(1/(10^pow)))) 2 (10^(pow+1)),Norm[residual[y]]},{y,yShadow-(1/(10^pow)),yShadow+(1/(10^pow)),(1/(2 (10^(pow+1))))}];

ListLinePlot[data,PlotStyle->Thickness[0.008], Mesh->All,MeshStyle->{Red,PointSize[0.015]},BaseStyle->14, Frame->True, PlotRange->{{0,Length[data]},{-.1, 3}} ,FrameTicks->{{{20,"yShadow"},{0,"yShad-(10^(-17))"},{40,"yShad+(10^(-17))"}}, Automatic, {},{}}]

Code

n=20;yShadow=seed[150];

pow=15;
data=Table[{-(y-(yShadow+(1/(10^pow)))) 20 (10^(pow+1)),Norm[residual[y]]},{y,yShadow-(1/(10^pow)),yShadow+(1/(10^pow)),(1/(20 (10^(pow+1))))}];

ListLinePlot[data,PlotStyle->Thickness[0.008], Mesh->All,MeshStyle->{Red,PointSize[0.01]},BaseStyle->14, Frame->True, PlotRange->{{0,Length[data]},{-.1, 6}} ,FrameTicks->{{{200,"yShadow"},{0,"yShad-(10^(-18))"},{400,"yShad+(10^(-17))"}}, Automatic, {},{}}]

Surprise: It appears that this works just as well without the Jacobian info. The Jacobian info is complicated.

Rob Knapp: Optimization

Abstract

Many dynamical systems are known to admit shadowing. This means that a noisy trajectory computed, say, using machine precision, is in fact the exact trajectory from a slightly different set of initial conditions. We will show how, in three very different examples, Mathematica's FindMinimum command can be used to minimize the least squares error and so find the exact values of the shadowing set of initial conditions.

The first case is the classic example of the quadratic map, (Subscript[f, r])=r x (1-x), where r is a value for which iterations are sensitively dependent on initial conditions, such as 4 or 38/10. In the latter case, for example, the following graph shows 100 iterations of (Subscript[f, 3.8]) using the starting value (Subscript[x, 0])=0.1 and machine precision. If one starts with
(Subscript[x, 0])= 0.100000000000000021043239304570322654494662756467086434194016030 and uses high precision, the resulting trajectory is within (10^(-15)) of the noisy one.

f[x_] := 3.8 x (1-x);
ListLinePlot[NestList[f, 0.1, 100], Frame->True, Axes->False, AspectRatio-> 0.4,FrameTicks->{Automatic, Automatic,None,None}, PlotStyle->Thickness[0.004], ImageSize->400]

A second example arises from Problem 2 of the SIAM 100-Digit Challenge, where a photon of light starts at ((1/2),(1/10)) and heads east at unit speed reflecting off circular mirrors of radius (1/3) centered at integer lattice points. Again, we can compute a machine precision path for the photon, and then find a high-precision value (Subscript[y, 0]) so that the noisy path (up to time 10) is the exact path for a photon starting at ((1/2),(1/10)+(Subscript[y, 0])) where we can specify (Subscript[y, 0]) exactly: 0.1000000000000000414132495877039710373799… .

In these two cases, getting derivative information for the iterated function is a key step, where the derivative is respect to change in the initial conditions, (Subscript[x, 0]) in the first case and (Subscript[y, 0]) in the second. A third example arises from the Duffing differential equation (x''-a x + b (x^3)+c x'=d cos t), where one can use the equation to discover the derivative information that one needs for shadowing.

The following image shows a typical machine-precision orbit (starts at (-1.,1.)) in the phase plane for the Duffing equation x''-x +  (x^3)+(15/100)x'=(3/10) cos t. The optimization method can find a nearby initial condition whose high-precision orbit matches the noisy one essentially exactly.

sol=x[t]/.First[NDSolve[{{x''[t]+(15/100) x'[t]-x[t]+((x[t])^3)==(3/10) Cos[t]},x[0]==-1,x'[0]==1},x,{t,100}]];
ParametricPlot[Evaluate[{sol, D[sol,t]}],{t,0,100},PlotStyle->Thick,Frame->True, Axes->False, FrameTicks->{Automatic, Automatic, None, None}]

Rough notes on why this works, by SW   Probably omit

I let x(t,α) denote the solution to the eqn starting from α.

dx0(t)=((x(t,(Subscript[x, 0]))-x(t,0))/(Subscript[x, 0]))   sensitivity ratio for starting at (Subscript[x, 0]), as (Subscript[x, 0])-> 0
Initial value:  dx0(0)=((x(0,(Subscript[x, 0]))-x(0,0))/(Subscript[x, 0]))=(((Subscript[x, 0])-0)/(Subscript[x, 0]))=1

General:  dx0''(t)=((x''(t,(Subscript[x, 0]))-x''(t,0))/(Subscript[x, 0]))=(((((-15)/100)x'(t,(Subscript[x, 0]))-x(t,(Subscript[x, 0]))+((x(t,(Subscript[x, 0])))^3)+(3/10)cos t) - (((-15)/100)x'(t,0)-x(t,0)+((x(t,0))^3))+(3/10)cos t))/(Subscript[x, 0]))=

 ((((-15)/100)(x'(t,(Subscript[x, 0]))-x'(t,0))-x(t,(Subscript[x, 0]))+((x(t,(Subscript[x, 0])))^3))+x(t,0)+((x(t,0))^3)))/(Subscript[x, 0]))=
 
 ((((-15)/100)(x'(t,(Subscript[x, 0]))-x'(t,0)))/(Subscript[x, 0]))+((-x(t,(Subscript[x, 0]))+x(t,0)+((x(t,(Subscript[x, 0])))^3)-((x(t,0))^3)))/(Subscript[x, 0]))= ((((-15)/100)(x'(t,(Subscript[x, 0]))-x'(t,0)))/(Subscript[x, 0]))-((x(t,(Subscript[x, 0]))-x(t,0))/(Subscript[x, 0]))+((x(((t,(Subscript[x, 0])))^3)-((x(t,0))^3))/(Subscript[x, 0]))=
 
  ((((-15)/100)(x'(t,(Subscript[x, 0]))-x'(t,0)))/(Subscript[x, 0]))+dx0(t)-((((x(t,(Subscript[x, 0])))^3)-((x(t,0))^3))/(Subscript[x, 0]))=

Now this last term is the sensitivity of (x^3). So standard calculus:   ((∂((x^3)))/(∂(Subscript[x, 0])))=3(x^2)((∂(x))/(∂(Subscript[x, 0])))

This notation is a little confusing. For this last step: Imagine we have f(t,(Subscript[x, 0])). dx0 represents how f(t,(Subscript[x, 0])) changes when (Subscript[x, 0]) changes. It is ((∂f)/(∂(Subscript[x, 0]))). But now suppose we have ((f(t,(Subscript[x, 0])))^3). Then we want ((∂((f(t, (Subscript[x, 0])))^3))/(∂(Subscript[x, 0]))), which is 3 ((f(t, (Subscript[x, 0])))^2) ((∂f)/(∂(Subscript[x, 0]))). That is what is going on here, where x(t,(Subscript[x, 0])) is used for f.

((((-15)/100)(x'(t,(Subscript[x, 0]))-x'(t,0)))/(Subscript[x, 0]))+dx0(t)-3((x(t,0))^2) dxo(t)=
        dx0''(t)=((-15)/100)dx0'(t)+dx0(t)-3((x(t))^2) dx0(t)

A Quick Experiment as a check, SW   Probably omit

Here we compute two solutions out to 100. Note: We have been looking at initial values (-1,1). Here we are switching to (0,0).

time=2;
sol1=x[t] /. First[NDSolve[Join[duff,{ x[0]==0, x'[0]==0}], x[t], {t,0,time}]];
sol1Close=x[t] /. First[NDSolve[Join[duff,{ x[0]==0.01, x'[0]==0}], x[t], {t,0,time}]];
{sol1/. t-> time, sol1Close /. t-> time}

%[[2]]-%[[1]]

So the rate is 3.19

threeSols={x[t],dx0[t],dxp0[t]}/.First[NDSolve[{duff,x[0]==0,x'[0]==0,sx0,sxp0},{x[t],dx0[t],dxp0[t]},{t,time}]]

threeSols[[2]]/.t->time

This is close to the observed value of 3.19 for the error-scale factor.

Here is a graph of the sensitivity function, log10s:

Plot[Log[10,Abs[threeSols[[2]]]],{t,0,time},PlotStyle->Thick]

What do we learn here:

The sensitivity function predicts how sensitive the soln at 2 is to a change in the init conditions. A factor of 3 shows up. As the 0.001 decreases the prediction is even better. And as we go out to 100 or so the sensitivity becomes very large: (10^8).

But we can compute this sensitivity accurately and with high precision. And we can compute it far out.

Attempt to do it same way as quadratic

ShadowDuffing[noisy_, n_, intivals_] := (
r[{x0_?NumericQ,xp0_?NumericQ}] :=Block[{$MinPrecision = n,$MaxPrecision = n},
Transpose[{x[samples],x'[samples]}/.hpssol[{x0,xp0}]]-Take[noisy,n+1]];

FindMinimum[0,{{x0, SetPrecision[xStart,n]},{xp0,SetPrecision[ xpStart,n]}}, WorkingPrecision->wp, AccuracyGoal->0.8 ag, PrecisionGoal->0.8 pg, 
Method->{"LevenbergMarquardt","Residual"->r[{x0,xp0}],"Jacobian"->J[{x0, xp0},n]}])

Do[sol=Shadow[noisy,nn,If[nn==50,1/10,x0/.sol[[2]]]],{nn,50,100,50}]
shadowSeed100=x0/.sol[[2]]

Photon Trajectory code

Clear[h];
norm[a_] :=(Sqrt[a.a]);

myRealDigits[a_Real] := Module[{man,exp,ss},
{man,exp}=MantissaExponent[Abs[a]];
man=ToString[InputForm[man]];
ss=StringPosition[man, "`"];
If[ss≠{}, man = StringDrop[man, {ss[[1,1]], StringLength[man]}]];
{ToExpression/@PadRight[Characters[StringDrop[man,2]],Ceiling[ Precision[a]]], exp}];

mySetPrecision[0|0.,_]:=0;
mySetPrecision[a_,prec_]:=(
rd=myRealDigits[SetPrecision[a,prec+1]];
fd=Take[rd[[1]],prec];
If[rd[[1,prec+1]]≥5,fd=ReplacePart[fd,fd[[-1]]+1,prec]];
ffd=FromDigits[{fd,rd[[2]]}];
((Sign[a] N[Numerator[ffd],prec])/(N[Denominator[ffd],prec])));

(* not necessary in version 5 *)
Unprotect[Round]; 
Round[Interval[{a_, b_}]] := Interval[{Round[a], Round[b]}];
Protect[Round]; 

iMin[{Interval[{a_, b_}], Interval[{c_, d_}]}]  :=  Interval[{Min[a,c], Min[b,d]}];

IntervalToReal[Interval[{a_,b_}]]:=((a+b)/2);
IntervalToReal[a_]:=a;
RealToInterval[r_,d_]:=Interval[r+{-1,1} d];

diam[Interval[{a_, b_}]] := b-a;
diam[{i__Interval}] := Max[diam /@ {i}];
h[{a_,b_}]:=9{{(b^2)-(a^2),-2 a b},{-2 a b,(a^2)-(b^2)}};

PhotonTrajectory::usage = "PhotonTrajectory[p, v, targettime (opts)] gives the set of reflection points and the first and last point. If interval arithmetic is used, the results are all intervals. I also gives points at one second intervals between bounces.";

StartIntervalPrecision::usage="StartIntervalPrecision is an option to PhotonTrajectory that sets the starting size of the intervals. A setting of automatic means that the accuracy goal is used for the initial size.";
IntervalPrecisionStep::usage="IntervalPrecisionStep is an option to PhotonTrajectory that sets the number of powers of 10 by which the interval size is decreased after a failure.";
IntervalArithmetic::usage="IntervalArithmetic is a possible choice to the Method option of PhotonTrajectory.";
FixedPrecision::usage="FixedPrecision is a possible choice to the Method option of PhotonTrajectory.";
SignificanceArithmetic::usage = "SignificanceArithmetic is a possible choice to the Method option of PhotonTrajectory. When using this method it is essential that the initial conditions be approximate real numbers.";
FixedWorkingPrecision::usage = "FixedWorkingPrecision is an option to PhotonTrajectory that sets the fixed precision to be used when the FixedPrecision method us chosen.";

Options[PhotonTrajectory] := {StartIntervalPrecision ->Automatic,AccuracyGoal -> 12, IntervalPrecisionStep -> 1, FixedWorkingPrecision->16, Trace ->False, Method -> SignificanceArithmetic};
 
PhotonTrajectory[p_,v_,targettime_, opts___Rule]:=Module[{ans,mat,pint,t,a,b},
{ag,stint, ips,met, fp, trQ} ={AccuracyGoal, StartIntervalPrecision, IntervalPrecisionStep,Method,  FixedWorkingPrecision, Trace} /. {opts} /. Options[PhotonTrajectory];

If[stint===Automatic, stint = ag];

Which[met === FixedPrecision, 

sp[a_] := a/. z_Real :> mySetPrecision[z,fp];
p1=sp[N[p,fp]];v1=sp[N[v,fp]];pdata={sp[p1]};
resttime=targettime;
While[resttime>0,m=Round[p1+((2 v1)/3)];
s=t/.Solve[norm[p1+t v1-m]^2==1/9,t];
t0=Min[Cases[s,_?Positive]];
If[t0<resttime,
p1=sp[p1+t0 v1];v1=sp[h[p1-m].v1],  
  t0=Min[resttime,(2/3)];p1=sp[p1+t0 v1]];
resttime=sp[resttime-t0];
   AppendTo[pdata,p1]];
 pdata,


met === IntervalArithmetic,
lastpoint=Interval[{-∞,∞}];intsize=stint;wp=Max[17,2+intsize];
While[diam[lastpoint]>(10^(-ag)),
pint=N[(RealToInterval[#,(10^(-intsize))]&)/@p,wp];
vint=N[(RealToInterval[#,(10^(-intsize))]&)/@v,wp];
pathdata={N[pint,wp]};
resttime=Interval[{targettime,targettime}];
While[resttime>0,m=Round[pint+((2 vint)/3)];
S=t/.(Solve[(pint+t vint-m).(pint+t vint-m)==1/9,t]);
If[FreeQ[S,Power[Interval[{_?Negative, _?Positive}],_]],
           T=Min[Cases[S,_?Positive]],Break[]];
Which[
T≤resttime,
  pint+=T vint;vint=h[pint-m].vint; resttime-=T,
T>resttime &&  resttime≥2/3,pint+=2/3 vint; resttime-=2/3,
T>resttime && resttime<2/3,pint+=resttime vint;resttime=0,
True, Break[]];
AppendTo[pathdata,pint];
If[Precision [ {resttime, pint, vint, T}]< ag, Break[]]];
intsize+=ips; wp =Max[17,intsize+2];
lastpoint=pint+Table[RealToInterval[-Max[Abs[resttime]],Max[Abs[resttime]]],{2}]];
If[trQ, Print[StringForm["Initial condition interval radius is (10^(-``)).",intsize]]];
pathdata,

met === SignificanceArithmetic,
p1=p;v1=v;
    pathdata={p1};resttime=targettime;
     While[resttime>0,c=Round[p1+(2/3)v1];time=Min[Cases[t/.Solve[norm[p1+t v1-c]^2==1/9,t],_?Positive]];
If[time<resttime,
   p1+=time v1;v1=h[p1-c].v1,
time=Min[resttime,2/3];p1+=time v1];
resttime-=time;
  AppendTo[pathdata,p1]];
pathdata]];

Junk

yShadow=(1/10)+(414132495877/(10^28));
ShadowTraj=PhotonTrajectory[N[{(1/2),yShadow},300],{1,0},41,WorkingPrecision->200];
n=50;
residual=Norm[Take[ShadowTraj,n]-Take[base, n]]

Graphics[{{
Pink,Thickness[0.013],Line[Take[ShadowTraj,n]]},
{mirrors,Thickness[0.004],Line[Take[base,n]], {Red, PointSize[Large], Point[base[[{1}]]]}}},PlotRange->{{-4.,1.},{-1.8,1.8}}, ImageSize->450]

Extra Code

Shadow[noisy_, n_, start_] := (
residual[x0_?NumericQ] :=Block[{$MinPrecision = n,$MaxPrecision = n},
fTraj[x0,n]-Take[noisy,n+1]];
FindMinimum[0,{x0,SetPrecision[start, n]}, WorkingPrecision->n,
Method->{"LevenbergMarquardt","Residual"->residual[x0],"Jacobian"->sensitivity[x0,n]}])

ShadowSeed[s_, nmax_] := (noisy=SetPrecision[fTraj[N[s],nmax],nmax];Do[sol=Shadow[noisy,nn,If[nn==50,s,x0/.sol[[2]]]],{nn,50,nmax,50}];
shadowSeed100=x0/.sol[[2]]);

ratsw[x_Real]:=(rd=RealDigits[x];rd[[1]].(10^(Range[-1,-16,-1])) (10^(rd[[2]])))

x0=1. (10^(-20));
ss=ShadowSeed[x0,100]

When getting the difference use relative error, and also get difference by rationalizing the noisy start, so as to force zeroes.

(((ss-ratsw[x0]))/x0)