John Starrett

178 Reputation

8 Badges

20 years, 22 days

MaplePrimes Activity


These are replies submitted by John Starrett

@Markiyan Hirnyk 

1. Yes, that was an unusual problem, and it turns out that when I increase the number of digits to 30, as Preben suggested, it fixes the out of order problem in stats fit, and makes Statistics Fit give a better polynomial fit , in terms of the goal of the whole project too.

2. Yes, you and preben showed that the out of order problem was fixed by your respective methods, and I am grateful for your assistance.

3. The "new" topic is the project where the out of order problem showed up. I am not using only polynomial basis functions, but many other basis functions to obtain differential equations with certain periodic orbit structures, as well as other properties. Even though there is no obvious reason to believe that low degree polynomials should be able to produce the paths on a Mobius band, degree 2 and 3 polynomials can indeed do it, which is interesting. 

Thank you for your help and comments.

@Preben Alsholm 

Yes, increasing the digits does help to reduce the error. In fact, it fixes the order error in stats fit also. Thank you for your attention and the suggestions.

@Markiyan Hirnyk 

Markiyan, I may not have as much time as you do. I am not ignoring this discussion, but between teaching, committees, grant meetings and so on, I have not had time to get back to Maple primes until now. I will try to upload the worksheet now, but if it fails, I may have to paset it in as I did before.

 

Just so you understand my criteria, I am interested in getting a fit that forces a differential equation to have a solution as close as possible to the original curves that generated the sampled data. Of course, there are many possible vector fields that can realise a set of a few curves as solutions, and they can vary greatly away from the particular solutions I am seeking. But when a least squares solution to my sampled data does not produce the original curves when given initial conditions on those curves, then either the functions fit to the data are inadequate, and I am asking too much of them, or the solution method is in error.

Now, it can also be that a solution method gives a fit with a larger error than is optimal, but produces what looks like a good solution, while a highly accurate fit to the data gives a diverging solution to the differential equations. In this case, I am asking too much of the functions, as would be the case if I asked for a system of ODEs with first degree polynomial right hand sides to reproduce the edges of  Mobius band.

So in that case, even though the less than optimal solution to the least squares problem produced a good differential equation, one that produced the right solution, whie a highly accurate method produced a diverging solution, the less good solution is preferred, because it gives me a DE that reproduces my data.

Thanks for all your attention. If my upload does not work, I will chack back later tonight or tomorrow and paste in the text box.MobiusFitMethodsComparison.mw   

@Preben Alsholm 

Thanks, Preben. The direct method you suggested gives good results. Here are two plots, for a degree three polynomial vector field. The original data is reproduced almost exactly, and the order does not matter.

@Markiyan Hirnyk 

Thank you for directing me to the package DirectSearch. Unfortunately, the results applied to the particular problem I outlined above give a badly diverging vector field, whereas stats fit gives good results.  

I have attached pictures of the output, which should be the center and edges of a Mobius band. The pictures on the left are from DirectSearch, and on the right using stats fit.

I will try version 2 of DirectFit.

John Starrett

@Markiyan Hirnyk 

I don't understand your "Is that all?" post. I sent an attachment that provided context for the problem, and am now adding your and Preban's suggestions. I will try attaching the file again. I will install DirectSearch and DirectSearch2 after I get back and see how they work. Was there an attached file? I cannot tell by looking

This worksheet builds sets of closed curves m[i](x(t), y(t), z(t)) in three space, finds their tangents, then fits polynomial functions f(x,y,z), g(x,y,z), h(x,y,z) to the x, y and z components of the tangent vectors to those curves, globalizing the local vector fields. These polynomials become the right hand sides of a system of three ODEs whose solutions should include the original curves that generated the vector fields.

"x' = "f(x,y,z)

y' = g(x,y,z)
z' = g(x,y,z)

 

In this particular worksheet, we are finding differential equations whose solutions should include paths along a Mobius band wrapping around the z axis.

restart:with(plots):with(DEtools):with(LinearAlgebra):with(Statistics):with(stats):with(statplots):with(linalg):

This matrix twists the band

R1:=<<cos(Pi*t),0,sin(Pi*t)>|<0,1,0>|<-sin(Pi*t),0,cos(Pi*t)>>:

This one wraps the band around the z axis.

R2:=<<cos(2*Pi*t),sin(2*Pi*t),0>|<-sin(2*Pi*t),cos(2*Pi*t),0|<0,0,1>>>:

V0 is the offset vector, determining the center of the band. The parameters nmax and width determine how many equally spaced orbits around the band are computed, and how wide the band is. The commented subscripted Vs are to determine if order makes a difference for any of the solution methods.

V0:=<4,0,0>:
nmax:=2:
width:=3:

for i from 0 to nmax do:
 V[i]:=<width/2-i*width/nmax,0,0>:
od:

V[0]:=<1,0,0>:
V[1]:=<-1,0,0>:
V[2]:=<0,0,0>:

This loop constructs each curve and its derivative.

for j from 0 to nmax do:
 Mob[j]:=(R2.(V0+R1.V[j])):
 dMob[j]:=<diff(Mob[j][1],t),diff(Mob[j][2],t),diff(Mob[j][3],t)>:
od:

This loop makes plots of the curves and their derivatives

for k from 0 to nmax do:
 P[k]:=spacecurve([Mob[k][1],Mob[k][2],Mob[k][3],t=0..1],color=black,numpoints=500,scaling=constrained):
 dP[k]:=spacecurve([dMob[k][1],dMob[k][2],dMob[k][3],t=0..1],color=black,numpoints=500,scaling=constrained):
od:
display(seq(P[n],n=0..nmax));
display(seq(dP[n],n=0..nmax));

This loop computes a sequence of sample points around each curve and each derivative curve, with kmin  and  kmax determining the number of points along the curves. This the data we will fit.

kmax:=15:
kmin:=0:
for m from 0 to nmax do;
Lx[m]:=evalf(seq(subs(t=k/kmax,Mob[m][1]),k=kmin..kmax)):
Ly[m]:=evalf(seq(subs(t=k/kmax,Mob[m][2]),k=kmin..kmax)):
Lz[m]:=evalf(seq(subs(t=k/kmax,Mob[m][3]),k=kmin..kmax)):
dLx[m]:=evalf(seq(subs(t=k/kmax,dMob[m][1]),k=kmin..kmax)):
dLy[m]:=evalf(seq(subs(t=k/kmax,dMob[m][2]),k=kmin..kmax)):
dLz[m]:=evalf(seq(subs(t=k/kmax,dMob[m][3]),k=kmin..kmax)):
od:

Plot the curves and the sample points. 

SP:=ScatterPlot3D(Matrix([[seq(Lx[n],n=0..nmax)],[seq(Ly[n],n=0..nmax)],[seq(Lz[n],n=0..nmax)]]),symbolsize=22);
display(seq(P[q],q=0..nmax),SP);
dSP:=ScatterPlot3D(Matrix([[seq(dLx[n],n=0..nmax)],[seq(dLy[n],n=0..nmax)],[seq(dLz[n],n=0..nmax)]]),symbolsize=22);
display(seq(dP[q],q=0..nmax),dSP);

#dSP:=ScatterPlot3D(Matrix([[seq(dLx[n],n=0..nmax)],[seq(dLy[n],n=0..nmax)],[seq(dLz[n],n=0..nmax)]]),symbolsize=22);
#display(seq(dP[q],q=0..nmax),dSP);

Here are the concatenated sequences of  data.

XL:=[seq(Lx[p],p=0..nmax)]:
YL:=[seq(Ly[p],p=0..nmax)]:
ZL:=[seq(Lz[p],p=0..nmax)]:
dXL:=[seq(dLx[p],p=0..nmax)]:
dYL:=[seq(dLy[p],p=0..nmax)]:
dZL:=[seq(dLz[p],p=0..nmax)]:
nops(XL);


Set up the polynomial to be fit.

PolyDegree:=3;
Poly:=`randpoly/monomials/dense`([x,y,z],PolyDegree);
Poly1:=add(cat(a,i)*Poly[i],i=1..nops(Poly));
Poly1Coeff:=cat(a,1..nops(Poly));

The old way using stats and fit, and the current method using Statistics and Fit

W1a:=fit[leastsquare[[x,y,z,w], w=Poly1, {Poly1Coeff}]]([XL,YL,ZL,dXL]);
W2a:=fit[leastsquare[[x,y,z,w], w=Poly1, {Poly1Coeff}]]([XL,YL,ZL,dYL]);
W3a:=fit[leastsquare[[x,y,z,w], w=Poly1, {Poly1Coeff}]]([XL,YL,ZL,dZL]);
W1:=Fit(Poly1,Transpose(Matrix([XL,YL,ZL])),Vector(dXL),[x,y,z],output=leastsquaresfunction);
W2:=Fit(Poly1,Transpose(Matrix([XL,YL,ZL])),Vector(dYL),[x,y,z],output=leastsquaresfunction);
W3:=Fit(Poly1,Transpose(Matrix([XL,YL,ZL])),Vector(dZL),[x,y,z],output=leastsquaresfunction);

Plot the surfaces to compare the two methods easily. Note that there is alomst no difference for the polyniomial fits of degree 1 and 2, but a big difference for 3

xrange:=3:
yrange:=3:
zrange:=3:
PW1a:=implicitplot3d(rhs(W1a)=0,x=-xrange..xrange,y=-yrange..yrange,z=-zrange..zrange,style=patchcontour):
PW2a:=implicitplot3d(rhs(W2a)=0,x=-xrange..xrange,y=-yrange..yrange,z=-zrange..zrange,style=patchcontour):
PW3a:=implicitplot3d(rhs(W3a)=0,x=-xrange..xrange,y=-yrange..yrange,z=-zrange..zrange,style=patchcontour):
PW1:=implicitplot3d(W1=0,x=-xrange..xrange,y=-yrange..yrange,z=-zrange..zrange,style=patchcontour):
PW2:=implicitplot3d(W2=0,x=-xrange..xrange,y=-yrange..yrange,z=-zrange..zrange,style=patchcontour):
PW3:=implicitplot3d(W3=0,x=-xrange..xrange,y=-yrange..yrange,z=-zrange..zrange,style=patchcontour):
display(Array([PW1a,PW2a,PW3a]));
display(Array([PW1,PW2,PW3]));

Two different versions of differential equations, and some initial conditions

dXa:=diff(x(t),t)=subs({x=x(t),y=y(t),z=z(t)},rhs(W1a)):
dYa:=diff(y(t),t)=subs({x=x(t),y=y(t),z=z(t)},rhs(W2a)):
dZa:=diff(z(t),t)=subs({x=x(t),y=y(t),z=z(t)},rhs(W3a)):
dX:=diff(x(t),t)=subs({x=x(t),y=y(t),z=z(t)},W1):
dY:=diff(y(t),t)=subs({x=x(t),y=y(t),z=z(t)},W2):
dZ:=diff(z(t),t)=subs({x=x(t),y=y(t),z=z(t)},W3):
r:='r':
for r from 0 to nmax do:
 IC[r]:=[x(0)=V0[1]+V[r][1],y(0)=V0[2]+V[r][2],z(0)=V0[3]+V[r][3]]:
od:
IC:=seq(IC[h],h=0..nmax);

Plot the two different solutions side by side. The equations ending in 'a' are computed using stats fit

DEP1a:=DEplot3d([dXa,dYa,dZa],[x(t),y(t),z(t)],t=0..1,[IC],linecolor=black, numpoints=3500, thickness=1,scaling=constrained):
DEP1:=DEplot3d([dX,dY,dZ],[x(t),y(t),z(t)],t=0..1,[IC],linecolor=black, numpoints=3500, thickness=1,scaling=constrained):

display(Array([DEP1a,DEP1]),view=[-5.01..6.01,-6.01..6.01,-2.01..2.01],orientation=[70,30]);

 


Download MobiusFitMethodsComparison.mw

.

@Markiyan Hirnyk 

This only compares stats fir and Stats Fit. I will add the methods outlined above after my meeting.

@Markiyan Hirnyk 

I've been working. The first weeks of class are quite hectic. I'll post an example worksheet with comparison calculations sometime today.

@Preben Alsholm 

I will be testing these with my data later today. Thanks very much for your work and feedback

John Starrett

@Markiyan Hirnyk 

I will be testing these with my data later today. Thanks very much for your work and feedback

John Starrett

@Carl Love 

 

Thanks Carl. I'll try to make a simplified worksheet that exhibits the bug.

@Carl Love 

 

I am taking collections of closed three dimensional curves, computing their tangents, then sampling the curves and their tangents to fit polynomials as the right hand side of differential equations. When the differential equations have identical (almost) solutions to the original curves, the fit is deemed a success. 

Unfortunately, while the Statistics Fit command does not seem to be sensitive to rearrangement of the order of sampled curves like stats fit was, it does not always produce good results. For instance, while the stats fit would return third degree polynomial functions that reproduced the original curves as solutions to differential equations, Statistics Fit failed badly by giving highly divergent solutions from identical initial conditions.

Surprisingly, it returns almost identical 1st 2nd and 4th degree polynomial fits that give good differential equation solutions. 

Thank you for your response, Carl, and to Preben also.

 

John Starrett

@Preben Alsho

I intended to attach it, but couldn't figure out how to do that. I didn't want to paste it in line, but I guess it couldn't hurt. In any case, since it is depracated, it will probably do no good to complain. 

Test of out of order bug for polynomial curve fitting. Putting the data in different orders results in different fits.
restart:with(plots):with(DEtools):with(LinearAlgebra):with(Statistics):with(stats):with(statplots):with(linalg):
Three components of x data
X1:= 4., 3.196471867, 1.177050981, -1.108686614, -2.677050983, -3.000000000, -2.177050982, -.7454153486, .6770509837, 1.657630100, 2.;
X2:=2., 1.657630099, .6770509820, -.7454153508, -2.177050984, -3.000000000, -2.677050982, -1.108686612, 1.177050984, 3.196471868, 4.;
X3:=3., 2.427050983, .9270509814, -.9270509826, -2.427050984, -3., -2.427050982, -.9270509802, .9270509838, 2.427050984, 3.;
x data arranged two different ways
XA:=[X1,X2,X3];
XB:=[X1,X3,X2];
Three components of y data
Y1:=0., 2.322372752, 3.622590434, 3.412186543, 1.944991388, -1.230620284*10^(-9), -1.581720127, -2.294152555, -2.083748665, -1.204338762, 1.640827046*10^(-9);
Y2:=0., 1.204338763, 2.083748665, 2.294152555, 1.581720125, -1.230620284*10^(-9), -1.944991391, -3.412186545, -3.622590432, -2.322372750, 3.281654092*10^(-9);
Y3:=0., 1.763355757, 2.853169550, 2.853169549, 1.763355757, -1.230620284*10^(-9), -1.763355759, -2.853169550, -2.853169549, -1.763355756, 2.461240569*10^(-9);
y data arranged two different ways
YA:=[Y1,Y2,Y3];
YB:=[Y1,Y3,Y2];
Three components of z data
Z1:=0., .3090169944, .5877852524, .8090169944, .9510565165, 1., .9510565163, .8090169941, .5877852522, .3090169936, -4.102067615*10^(-10);
Z2:=-0., -.3090169944, -.5877852524, -.8090169944, -.9510565165, -1., -.9510565163, -.8090169941, -.5877852522, -.3090169936, 4.102067615*10^(-10);
Z3:=0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.;
z data arranged two different ways
ZA:=[Z1,Z2,Z3];
ZB:=[Z1,Z3,Z2];
Three components values at (x, y, z)
W1:=-0., -15.37729651, -23.33203217, -20.65400220, -9.803534311, 3.141592662, 12.35544766, 15.19998378, 12.52195383, 6.781685450, -9.020917841*10^(-9);
W2:=0., -6.781685456, -12.52195383, -15.19998379, -12.35544765, -3.141592646, 9.803534333, 20.65400220, 23.33203215, 15.37729650, -2.190794333*10^(-8);
W3:=-0., -11.07949099, -17.92699300, -17.92699299, -11.07949098, 7.732215294*10^(-9), 11.07949100, 17.92699300, 17.92699299, 11.07949097, -1.546443058*10^(-8);
values at (x,y,z) arranged two different ways
WA:=[W1,W2,W3];
WB:=[W1,W3,W2];
Set up a polynomial of degree three
PolyDegree:=3;
Poly:=`randpoly/monomials/dense`([x,y,z],PolyDegree);
Poly1:=add(cat(a,i)*Poly[i],i=1..nops(Poly));
Poly1Coeff:=cat(a,1..nops(Poly));
Fit the two identical data sets, except for the switched order, obtaining two different polynomials with non-trivial differences
FA:=fit[leastsquare[[x,y,z,w], w=Poly1, {Poly1Coeff}]]([XA,YA,ZA,WA]);
FB:=fit[leastsquare[[x,y,z,w], w=Poly1, {Poly1Coeff}]]([XB,YB,ZB,WB]);

To make it more clear, we plot the solutions implicitly
PA:=implicitplot3d(rhs(FA)=1,x=-4..4,y=-4..4,z=-4..4,grid=[40,40,40],orientation=[0,90],style=patchnogrid);
PB:=implicitplot3d(rhs(FB)=1,x=-4..4,y=-4..4,z=-4..4,grid=[40,40,40],orientation=[0,90],style=patchnogrid);
P:=implicitplot3d((rhs(FB)-rhs(FA))=1,x=-4..4,y=-4..4,z=-4..4,grid=[40,40,40],orientation=[0,90],style=patchnogrid);
and display them
display(PA);
display(PB);
display(P);

That did the trick. Thanks for all the help. Here is an example

http://www.aurapiercing.com/LogisticAndStableManifolds.eps

John Starrett

Is there a trick for making this in color?

John Starrett

1 2 3 Page 1 of 3