Lesson 15: Strange attractor restart; with(plots): N1:= proc(X) local X1, t1, t5, t7, t12, t16, t22, t24; X1:= X[1]; t1 := cos(X[2]); t5 := X[2]^2; t7 := X[1]^2; t12 := 1/(2.*t1*X[1]*X[2]+t1*t5-t7-2.*X[1]*X[2]); t16 := t7*X[2]+X[1]*t5-2.; t22 := sin(X[2]); t24 := t12*(X[1]+t22-2.); X[1]:= X[1]-t1*t12*t16+X[1]*(X[1]+2.*X[2])*t24; X[2]:= X[2]+t12*t16-X[2]*(2.*X1+X[2])*t24; end proc: NC:= Compiler[Compile](N1);

Looking around at what happens with different starting points, I found large regions that were attracted to one or the other of our two solutions, but also other regions where more complicated things happened. Here's one. I took this starting point: W[0]:= Vector(<1.0, -0.1>, datatype=float[8]); Now I want to use NC, my compiled version of Newt. Remember, this operates by changing its input. In order to do this, and still save the old value, I'll have to operate on a copy of the value. Now something like this won't quite work: W[1]:= W[0]; NC(W[1]); W[0], W[1]; The problem was that assigning Vectors works a bit differently from ordinary expressions. What was assigned to W[1] was the Vector W[0], not just the value of W[0]. So when NC changed W[1], that also affected W[0]. Instead, we can use copy which will produce a copy of our Vector. W[0]:= Vector(<1.0, -0.1>, datatype=float[8]): W[1]:= copy(W[0]): NC(W[1]): W[0], W[1]; for count from 1 to 20 do W[count]:= copy(W[count-1]); NC(W[count]); end do: W[16],W[17],W[18],W[19],W[20]; The last few iterations look like they might be starting to settle down to a 2-cycle, but it doesn't happen. I'll give the process some more time to settle down, then save and plot the next 40000 points. I don't need to save the first 1000 points. X:= copy(W[20]): for count from 1 to 1000 do NC(X) end do: Now to plot the next 50000 points, using pointplot from the plots package. You can give pointplot a Matrix whose columns are the points, with datatype=float[8]. Note how I can assign a column of the Matrix to be the Vector X (I don't need to use copy here). M:= Matrix(2,50000, datatype=float[8]): for count from 1 to 50000 do NC(X); M[1..2,count]:= X; end do: pointplot(M, symbol=point); PP:= %: The points appear to lie on two curves that are folded over on themselves. On closer examination, it turns out that there seem to be an infinite number of "folds'': this is an example of a so-called "strange attractor''. display(PP,view=[.0085 .. .0095, 14.5 .. 15]); display([PP,plottools[rectangle]([.0085,14.5],[.0095,15],colour=red)]); display(PP,view=[.00898 .. .00901, 14.56 .. 14.59]); display([PP,plottools[rectangle]([.00898,14.56],[.00901,14.59],colour=red)], view=[.0085 .. .0095, 14.5 .. 15]); For a closer look at what happens here, I'll follow 1001 initial points near the end of one of the curves (say on the line [.010, 13.06] to [.01, 13.07]), and see what happens to them as we iterate our Newton function. This is using a different version of animation: here you produce each frame of the animation separately, and then give display a list of the frames, with the option insequence=true. Each of the 1001 points will have a slightly different colour, going through the rainbow from red to violet. Those can be specified as COLOUR(HUE,c) where c goes from 0 to 1 (but 1 is the same as 0) M:= Matrix(2,1001, datatype=float[8]): for ii from 1 to 1001 do M[1,ii]:= 0.01 ; M[2,ii]:= 13.06 + 0.00001*(ii-1); end do: colours:= evalf([seq(COLOUR(HUE,i/1100),i=0..1000)]): frame[0]:= pointplot(M, symbol=point, colour=colours): frame[0]; Now for each frame we apply NC to all the points, and plot the new positions of all the points. for count from 1 to 60 do for ii from 1 to 1001 do X:= M[1..2,ii]; NC(X); M[1..2,ii]:= X; end do; frame[count]:= pointplot(M, symbol=point,colour=colours); end do: display([seq(frame[count],count=0..60)], insequence=true); The points start out very close together, but spread out to fill half of the attractor (and alternate between the two halves). One of the hallmarks of chaos is sensitive dependence on initial conditions: a very small difference in the initial point produces large differences after many iterations.

restart; The main integration command in Maple is int. It can be used for both antiderivatives and definite integrals: int(x^2, x); int(x^2, x = 1 .. 4); Iterated integrals can be done: int(int(x^2*y, x = 0 .. y), y = 0 .. 1); There is also an inert integration command Int, for when you want to write down an integral formula but don't want to actually do the integration. Int(x^2, x = 1 .. 4); A formula containing Int can be manipulated in various ways. When you want to actually do the integration, you can use the value command on any expression containing the Int. % = value(%);

As you know, it can be hard to find an antiderivative. Maple is pretty good at it, but not perfect. Here's a hard one it can do. F := int(sqrt(tan(x)), x); simplify(F); Should we believe this? The way to check an antiderivative is to differentiate it. dF := diff(F, x); It doesn't look much like LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkmbXNxcnRHRiQ2Iy1GIzYmLUkjbWlHRiQ2JVEkdGFuRicvJSdpdGFsaWNHUSZmYWxzZUYnLyUsbWF0aHZhcmlhbnRHUSdub3JtYWxGJy1JI21vR0YkNi1RMCZBcHBseUZ1bmN0aW9uO0YnRjcvJSZmZW5jZUdGNi8lKnNlcGFyYXRvckdGNi8lKXN0cmV0Y2h5R0Y2LyUqc3ltbWV0cmljR0Y2LyUobGFyZ2VvcEdGNi8lLm1vdmFibGVsaW1pdHNHRjYvJSdhY2NlbnRHRjYvJSdsc3BhY2VHUSYwLjBlbUYnLyUncnNwYWNlR0ZOLUkobWZlbmNlZEdGJDYkLUYjNiQtRjE2JVEieEYnL0Y1USV0cnVlRicvRjhRJ2l0YWxpY0YnRjdGN0Y3Rjc=, but maybe it simplifies. simplify(dF ); To check if LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2JVEhRicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUYjNiZGKy1GIzYmLUYsNiVRIkFGJ0YvRjItSSNtb0dGJDYtUSI9RicvRjNRJ25vcm1hbEYnLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0ZELyUpc3RyZXRjaHlHRkQvJSpzeW1tZXRyaWNHRkQvJShsYXJnZW9wR0ZELyUubW92YWJsZWxpbWl0c0dGRC8lJ2FjY2VudEdGRC8lJ2xzcGFjZUdRLDAuMjc3Nzc3OGVtRicvJSdyc3BhY2VHRlMtRiw2JVEiQkYnRi9GMkZARitGQEYrRkA=, rather than hoping Maple can simplify A to B, a better idea might be to try to simplify LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2JVEhRicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUYjNiZGKy1GIzYmLUYsNiVRIkFGJ0YvRjItSSNtb0dGJDYtUSgmbWludXM7RicvRjNRJ25vcm1hbEYnLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0ZELyUpc3RyZXRjaHlHRkQvJSpzeW1tZXRyaWNHRkQvJShsYXJnZW9wR0ZELyUubW92YWJsZWxpbWl0c0dGRC8lJ2FjY2VudEdGRC8lJ2xzcGFjZUdRLDAuMjIyMjIyMmVtRicvJSdyc3BhY2VHRlMtRiw2JVEiQkYnRi9GMkZARitGQEYrRkA= to 0 (or LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2JVEhRicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUYjNiQtSSZtZnJhY0dGJDYoLUYjNiQtRiw2JVEiQUYnRi9GMi9GM1Enbm9ybWFsRictRiM2JC1GLDYlUSJCRidGL0YyRj8vJS5saW5ldGhpY2tuZXNzR1EiMUYnLyUrZGVub21hbGlnbkdRJ2NlbnRlckYnLyUpbnVtYWxpZ25HRksvJSliZXZlbGxlZEdRJmZhbHNlRidGP0YrRj8= to 1). simplify(dF - sqrt(tan(x))); This won't always work: checking whether a complicated expression is 0 is not an easy task. Another possibility is to plot the difference. plot(dF - sqrt(tan(x)), x = 0 .. Pi/2 - 0.0001, axes = box); The differences are small: presumably due to roundoff error. This one might be harder: J := Int(sqrt(sec(x)^2 + 4), x); value(J); It's not just complicated, but it involves some special functions (EllipticPi and EllipticF). Is it correct? simplify(diff(value(J),x)-sqrt(sec(x)^2 + 4)); But actually there's a relatively simple, elementary solution that Maple misses (I'm not sure exactly why). There is a Student[Calculus1] package that can help with integration at the Math 101 level. You can, for example, ask it for a hint on how to do an integral. with(Student[Calculus1]): Hint(J); It's telling us to try the change of variables LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2JVEhRicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUYjNictRiw2JVEidUYnRi9GMi1JI21vR0YkNi1RIj1GJy9GM1Enbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRkIvJSlzdHJldGNoeUdGQi8lKnN5bW1ldHJpY0dGQi8lKGxhcmdlb3BHRkIvJS5tb3ZhYmxlbGltaXRzR0ZCLyUnYWNjZW50R0ZCLyUnbHNwYWNlR1EsMC4yNzc3Nzc4ZW1GJy8lJ3JzcGFjZUdGUS1GIzYmLUYsNiVRJHRhbkYnL0YwRkJGPi1GOzYtUTAmQXBwbHlGdW5jdGlvbjtGJ0Y+RkBGQ0ZFRkdGSUZLRk0vRlBRJjAuMGVtRicvRlNGaG4tSShtZmVuY2VkR0YkNiQtRiM2JC1GLDYlUSJ4RidGL0YyRj5GPkY+RitGPkYrRj4=. It will also do that for you. J1:= Rule[%](J); What next? Hint(%); One of these looks interesting, but I'd prefer a different letter s instead of the name u1 Maple chose for the new variable. JSFH J2:=Rule[change,s=sqrt(5+u^2)-u,s](J1); OK, that's a rational function. We might tackle that with partial fractions. But let's just let Maple do it. value(rhs(%)); Now we have to substitute back to get the integral in terms of u, and then the original variable x. eval(%,s = sqrt(u^2+5)-u); G:= eval(%,u = tan(x)); That's our answer. Let's check it. dG:= diff(G,x); simplify(dG - sqrt(sec(x)^2+4)); Some antiderivatives just can't be expressed in "closed form" int(exp(sin(x)), x); Others need special functions. This one is important in probability. int(exp(-x^2/2), x);

copy insequence=true (option for display to produce an animation) COLOUR(HUE,...)) (option for colour in plot commands) int Int value indets lprint Student[Calculus1] package Hint (in Student[Calculus1]) Rule (in Student[Calculus1]) change (in Student[Calculus1]) JSFH

LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkjbWlHRiQ2JVEhRicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0Yn