Lesson 13: Newton in several variables restart;

JSFH There are very big and capable LinearAlgebra and VectorCalculus packages that can do lots of things with Vectors and Matrices, but we won't need to go much beyond the very basic level of Vector and Matrix operations from Lesson 12. I'll need one command from the VectorCalculus package: Jacobian. Consider a column Vector V whose components are expressions depending on several variables. V := <f(x,y,z), g(x,y,z), h(x,y,z)>; The Jacobian of F is the Matrix of partial derivatives of the components of F with respect to each of the variables. Since this is the only member of the VectorCalculus package I'll want to use today, and I want to avoid some side-effects of VectorCalculus, I'll only take this one procedure from VectorCalculus. You can do that with an extra input to with. with(VectorCalculus, Jacobian); J := Jacobian(V,[x,y,z]); Notice that the rows of J correspond to the components of V, and the columns correspond to the variables. You could think of each row of J as the gradient vector of a component of V. It will be convenient for us to consider functions from Vectors to Vectors. You can define these like this: F := X -> <X[1] + X[2]*X[3], X[2] - X[1]^2, cos(X[3])>; F(<a,b,c>); Unfortunately there isn't an equivalent of the D operator to differentiate such a function: Jacobian has to be applied to a Vector of expressions, not to a function such as F. What you can do is this: Jacobian(F(X),[X[1],X[2],X[3]]); And then you can make this into a function using unapply. J:= unapply(%,X); J(<a,b,c>); J(<1,2,3>); Why am I using unapply here rather than ->? badJ:= X -> Jacobian(F(X),[X[1],X[2],X[3]]); This would work if the components of X were variables: badJ(<x,y,z>); but not with constants: badJ(<1,2,3>); The Jacobian command needs to take an expression involving variables and differentiate it with respect to variables. Here badJ is calling Jacobian with a constant vector F(<1, 2, 3>) and a list of constants [1, 2, 3]. Maple complains, and rightly so. On the other hand, my definition of J computed the Jacobian once, using variables, and then used that to define a function that can be given either variables or constants.

So far we've seen Newton's method used for solving one equation in one variable. Now I want to look at the extension of this to solving a system of equations in several variables. The basic idea of Newton's method is the same: Start with an initial guess . Consider the linear approximation to LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2Ji1GLDYlUSJGRicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUkjbW9HRiQ2LVEwJkFwcGx5RnVuY3Rpb247RicvRjhRJ25vcm1hbEYnLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0ZCLyUpc3RyZXRjaHlHRkIvJSpzeW1tZXRyaWNHRkIvJShsYXJnZW9wR0ZCLyUubW92YWJsZWxpbWl0c0dGQi8lJ2FjY2VudEdGQi8lJ2xzcGFjZUdRJjAuMGVtRicvJSdyc3BhY2VHRlEtSShtZmVuY2VkR0YkNiQtRiM2JC1GLDYlUSJYRidGNEY3Rj5GPkY+RitGPg== near the initial guess. Find the point where the linear approximation is 0. Use that as the next guess. Repeat this process until it appears to converge. The linear approximation to LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2Ji1GLDYlUSJGRicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUkjbW9HRiQ2LVEwJkFwcGx5RnVuY3Rpb247RicvRjhRJ25vcm1hbEYnLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0ZCLyUpc3RyZXRjaHlHRkIvJSpzeW1tZXRyaWNHRkIvJShsYXJnZW9wR0ZCLyUubW92YWJsZWxpbWl0c0dGQi8lJ2FjY2VudEdGQi8lJ2xzcGFjZUdRJjAuMGVtRicvJSdyc3BhY2VHRlEtSShtZmVuY2VkR0YkNiQtRiM2JC1GLDYlUSJYRidGNEY3Rj5GPkY+RitGPg== near LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2Ji1GLDYlUSJYRicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUkjbW9HRiQ2LVEiPUYnL0Y4USdub3JtYWxGJy8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGQi8lKXN0cmV0Y2h5R0ZCLyUqc3ltbWV0cmljR0ZCLyUobGFyZ2VvcEdGQi8lLm1vdmFibGVsaW1pdHNHRkIvJSdhY2NlbnRHRkIvJSdsc3BhY2VHUSwwLjI3Nzc3NzhlbUYnLyUncnNwYWNlR0ZRLUYsNiVRIkFGJ0Y0RjdGPkYrRj4= is LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYqLUkjbWlHRiQ2JVEiRkYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JKG1mZW5jZWRHRiQ2JC1GIzYkLUYsNiVRIkFGJ0YvRjIvRjNRJ25vcm1hbEYnRj0tSSNtb0dGJDYtUSIrRidGPS8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGRS8lKXN0cmV0Y2h5R0ZFLyUqc3ltbWV0cmljR0ZFLyUobGFyZ2VvcEdGRS8lLm1vdmFibGVsaW1pdHNHRkUvJSdhY2NlbnRHRkUvJSdsc3BhY2VHUSwwLjIyMjIyMjJlbUYnLyUncnNwYWNlR0ZULUYsNiVRIkpGJ0YvRjJGNS1GQDYtUSIuRidGPUZDRkZGSEZKRkxGTkZQL0ZTUSYwLjBlbUYnL0ZWRmhuLUY2NiQtRiM2Ji1GLDYlUSJYRidGL0YyLUZANi1RKCZtaW51cztGJ0Y9RkNGRkZIRkpGTEZORlBGUkZVRjpGPUY9Rj0= where LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2Ji1GLDYlUSJKRicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUkjbW9HRiQ2LVEwJkFwcGx5RnVuY3Rpb247RicvRjhRJ25vcm1hbEYnLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0ZCLyUpc3RyZXRjaHlHRkIvJSpzeW1tZXRyaWNHRkIvJShsYXJnZW9wR0ZCLyUubW92YWJsZWxpbWl0c0dGQi8lJ2FjY2VudEdGQi8lJ2xzcGFjZUdRJjAuMGVtRicvJSdyc3BhY2VHRlEtSShtZmVuY2VkR0YkNiQtRiM2JC1GLDYlUSJBRidGNEY3Rj5GPkY+RitGPg== is the Jacobian matrix for LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2Ji1GLDYlUSJGRicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUkjbW9HRiQ2LVEwJkFwcGx5RnVuY3Rpb247RicvRjhRJ25vcm1hbEYnLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0ZCLyUpc3RyZXRjaHlHRkIvJSpzeW1tZXRyaWNHRkIvJShsYXJnZW9wR0ZCLyUubW92YWJsZWxpbWl0c0dGQi8lJ2FjY2VudEdGQi8lJ2xzcGFjZUdRJjAuMGVtRicvJSdyc3BhY2VHRlEtSShtZmVuY2VkR0YkNiQtRiM2JC1GLDYlUSJYRidGNEY3Rj5GPkY+RitGPg== at LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYnLUkjbWlHRiQ2JVEiWEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RIj1GJy9GM1Enbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRj0vJSlzdHJldGNoeUdGPS8lKnN5bW1ldHJpY0dGPS8lKGxhcmdlb3BHRj0vJS5tb3ZhYmxlbGltaXRzR0Y9LyUnYWNjZW50R0Y9LyUnbHNwYWNlR1EsMC4yNzc3Nzc4ZW1GJy8lJ3JzcGFjZUdGTC1GLDYlUSJBRidGL0YyLUY2Ni1RIi5GJ0Y5RjtGPkZARkJGREZGRkgvRktRJjAuMGVtRicvRk5GVkY5 Thus in two variables: F:= X -> <f(X[1],X[2]), g(X[1],X[2])>; J := unapply(Jacobian(F(X),[X[1],X[2]]),X); F(<a[1],a[2]>) + J(<a[1],a[2]>) . (<x[1],x[2]>-<a[1],a[2]>); Now to find where this is 0. Assuming the Matrix J(A) is invertible, the solution of LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYuLUkjbWlHRiQ2JVEiRkYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JKG1mZW5jZWRHRiQ2JC1GIzYkLUYsNiVRIkFGJ0YvRjIvRjNRJ25vcm1hbEYnRj0tSSNtb0dGJDYtUSIrRidGPS8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGRS8lKXN0cmV0Y2h5R0ZFLyUqc3ltbWV0cmljR0ZFLyUobGFyZ2VvcEdGRS8lLm1vdmFibGVsaW1pdHNHRkUvJSdhY2NlbnRHRkUvJSdsc3BhY2VHUSwwLjIyMjIyMjJlbUYnLyUncnNwYWNlR0ZULUYsNiVRIkpGJ0YvRjJGNS1GQDYtUTEmSW52aXNpYmxlVGltZXM7RidGPUZDRkZGSEZKRkxGTkZQL0ZTUSYwLjBlbUYnL0ZWRmhuLUY2NiQtRiM2Ji1GLDYlUSJYRidGL0YyLUZANi1RKCZtaW51cztGJ0Y9RkNGRkZIRkpGTEZORlBGUkZVRjpGPUY9LUZANi1RIn5GJ0Y9RkNGRkZIRkpGTEZORlBGZ25GaW4tRkA2LVEiPUYnRj1GQ0ZGRkhGSkZMRk5GUC9GU1EsMC4yNzc3Nzc4ZW1GJy9GVkZbcEZkby1JI21uR0YkNiRRIjBGJ0Y9Rj0= is LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYyLUkjbWlHRiQ2JVEiWEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RIn5GJy9GM1Enbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRj0vJSlzdHJldGNoeUdGPS8lKnN5bW1ldHJpY0dGPS8lKGxhcmdlb3BHRj0vJS5tb3ZhYmxlbGltaXRzR0Y9LyUnYWNjZW50R0Y9LyUnbHNwYWNlR1EmMC4wZW1GJy8lJ3JzcGFjZUdGTC1GNjYtUSI9RidGOUY7Rj5GQEZCRkRGRkZIL0ZLUSwwLjI3Nzc3NzhlbUYnL0ZORlNGNS1GLDYlUSJBRidGL0YyRjUtRjY2LVEqJnVtaW51czA7RidGOUY7Rj5GQEZCRkRGRkZIL0ZLUSwwLjIyMjIyMjJlbUYnL0ZORmZuRjUtRiw2JVEiSkYnRi9GMi1JJW1zdXBHRiQ2JS1JKG1mZW5jZWRHRiQ2JC1GIzYkRlVGOUY5LUYjNiVGWC1JI21uR0YkNiRRIjFGJ0Y5RjkvJTFzdXBlcnNjcmlwdHNoaWZ0R1EiMEYnRjUtRjY2LVEiLkYnRjlGO0Y+RkBGQkZERkZGSEZKRk1GNS1GLDYlUSJGRidGL0YyRl5vRjk= Let's try an example. sys := [x^2 + x*y*z + y^2*z^3 = 3, x*y^2 - y*z^2 - 2*x^2 = 0, x^2*y + y^2*z + z^4 = 4]; I want to write this as LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2JkYrLUYjNidGKy1GIzYmLUYsNiVRIkZGJy8lJ2l0YWxpY0dRJXRydWVGJy8lLG1hdGh2YXJpYW50R1EnaXRhbGljRictSSNtb0dGJDYtUTAmQXBwbHlGdW5jdGlvbjtGJy9GPFEnbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRkYvJSlzdHJldGNoeUdGRi8lKnN5bW1ldHJpY0dGRi8lKGxhcmdlb3BHRkYvJS5tb3ZhYmxlbGltaXRzR0ZGLyUnYWNjZW50R0ZGLyUnbHNwYWNlR1EmMC4wZW1GJy8lJ3JzcGFjZUdGVS1JKG1mZW5jZWRHRiQ2JC1GIzYkLUYsNiVRIlhGJ0Y4RjtGQkZCRkItRj82LVEiPUYnRkJGREZHRklGS0ZNRk9GUS9GVFEsMC4yNzc3Nzc4ZW1GJy9GV0Zeby1JI21uR0YkNiRRIjBGJ0ZCRkJGK0ZCRitGQg==, where F is a vector-valued function. I'll subtract the right sides from the left: map( lhs - rhs, sys); Change x, y, z to the components of the vector X. eval(%,{x=X[1],y=X[2],z=X[3]}); Make this into a column Vector. Vector(%); Use unapply to make this into a function of X. F:= unapply(%,X); Now to get the Jacobian: J := unapply(Jacobian(F(X),[X[1],X[2],X[3]]),X); The Newton iteration will be Newt:= X -> evalf(X - J(X)^(-1) . F(X)); For example, if we take this starting point W[0]:= <1, 2, 3>; Newt(W[0]); Did that bring us closer to a solution? F(%); F(W[0]); Maybe... Let's try a few more iterations. for count from 1 to 10 do W[count]:= Newt(W[count-1]) end do; It seems to have pretty much settled down (except for roundoff error). F(W[10]); So it seems there is a solution here. Is it the only one? Let's try a different initial condition. Y[0]:= <0,0,1>; for count from 1 to 10 do Y[count]:= Newt(Y[count-1]) end do; Oops... LinearAlgebra[Determinant](J(<x,y,z>)); Note each term has at least two of the variables x, y, z. So to get the Jacobian to be nonsingular, you need at least two of x, y, z to be nonzero. Try again... Y[0]:= <1,0,1>; for count from 1 to 10 do Y[count]:= Newt(Y[count-1]) end do; Let's try to get a picture. with(plots): The implicitplot3d command in the plots package is similar to implicitplot, but in three dimensions instead of two, and without gridrefine or crossingrefine options. It does have the grid option (with a list of three integers), but high numbers greatly increase the time needed to plot. implicitplot3d(sys,x=-10..10,y=-10..10,z=-10..10,colour=[red,blue,green],axes=box); implicitplot3d(sys,x=0..2,y=0..2,z=0..1,colour=[red,blue,green],axes=box,style=patchnogrid,lightmodel=light4); Here's my small contribution to Maple. The intersectplot command, also in the plots package, plots in 3 dimensions the curve where two surfaces intersect. Here's the intersection of the first two surfaces. intersectplot(sys[1],sys[2],x=-10..10,y=-10..10,z=-10..10,colour=red, axes=box); Here are the intersections of all three pairs of surfaces. At the intersection points, all three curves should meet. display([%,intersectplot(sys[2],sys[3],x=-10..10,y=-10..10,z=-10..10,colour=blue), intersectplot(sys[1],sys[3],x=-10..10,y=-10..10,z=-10..10,colour=green)],axes=box); IP1:= %: It would be nice to see the solutions we found earlier in here. The pointplot3d command in the plots package can be used to plot a list of points in 3 dimensions. While we're at it, I'll also zoom in a bit using the view option. display([IP1, pointplot3d([W[10],Y[10]],symbol=solidbox, symbolsize=10, colour=black)],view=[-3..3,-3..3,-3..3]); How many solutions are there in all? Newton won't tell us. Maybe solve will. solve(sys); Somewhat surprisingly, we need Digits to be pretty high to get reasonably accurate numerical values of the solutions. Even 20 wasn't enough. evalf([allvalues(%)],25); sols:=remove(has,%,I); pts:= [seq(eval(<x,y,z>,sols[i]), i = 1..nops(sols))]; W[10], Y[10]; display([IP1, pointplot3d(pts,symbol=solidbox, symbolsize=10, colour=black)],view=[-3..3,-3..3,-3..3]);

We've seen several ways of plotting curves in Maple. For curves in the LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEjeHlGJy8lJ2l0YWxpY0dRJXRydWVGJy8lLG1hdGh2YXJpYW50R1EnaXRhbGljRicvRjNRJ25vcm1hbEYn planeLUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkjbWlHRiQ2I1EhRic=of the form LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYpLUkjbWlHRiQ2JVEieUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RIn5GJy9GM1Enbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRj0vJSlzdHJldGNoeUdGPS8lKnN5bW1ldHJpY0dGPS8lKGxhcmdlb3BHRj0vJS5tb3ZhYmxlbGltaXRzR0Y9LyUnYWNjZW50R0Y9LyUnbHNwYWNlR1EmMC4wZW1GJy8lJ3JzcGFjZUdGTC1GNjYtUSI9RidGOUY7Rj5GQEZCRkRGRkZIL0ZLUSwwLjI3Nzc3NzhlbUYnL0ZORlNGNS1GLDYlUSJmRidGL0YyLUkobWZlbmNlZEdGJDYkLUYjNiQtRiw2JVEieEYnRi9GMkY5RjlGOQ== you can use plot(f(x), x=a..b); For curves given by an implicit equation such as LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2JkYrLUYjNidGKy1GIzYmLUYsNiVRImZGJy8lJ2l0YWxpY0dRJXRydWVGJy8lLG1hdGh2YXJpYW50R1EnaXRhbGljRictSSNtb0dGJDYtUTAmQXBwbHlGdW5jdGlvbjtGJy9GPFEnbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRkYvJSlzdHJldGNoeUdGRi8lKnN5bW1ldHJpY0dGRi8lKGxhcmdlb3BHRkYvJS5tb3ZhYmxlbGltaXRzR0ZGLyUnYWNjZW50R0ZGLyUnbHNwYWNlR1EmMC4wZW1GJy8lJ3JzcGFjZUdGVS1JKG1mZW5jZWRHRiQ2JC1GIzYmLUYsNiVRInhGJ0Y4RjstRj82LVEiLEYnRkJGRC9GSEY6RklGS0ZNRk9GUUZTL0ZXUSwwLjMzMzMzMzNlbUYnLUYsNiVRInlGJ0Y4RjtGQkZCRkItRj82LVEiPUYnRkJGREZHRklGS0ZNRk9GUS9GVFEsMC4yNzc3Nzc4ZW1GJy9GV0Znby1JI21uR0YkNiRRIjBGJ0ZCRkJGK0ZCRitGQg== you can use implicitplot in the plots package. implicitplot(f(x,y)=0, x=a..b, y=c..d); A third kind of curve is a parametric curve, where x and y are given as functions of a parameter t, e.g. LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2JkYrLUYjNihGKy1GIzYnLUYsNiVRInhGJy8lJ2l0YWxpY0dRJXRydWVGJy8lLG1hdGh2YXJpYW50R1EnaXRhbGljRictSSNtb0dGJDYtUSI9RicvRjxRJ25vcm1hbEYnLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0ZGLyUpc3RyZXRjaHlHRkYvJSpzeW1tZXRyaWNHRkYvJShsYXJnZW9wR0ZGLyUubW92YWJsZWxpbWl0c0dGRi8lJ2FjY2VudEdGRi8lJ2xzcGFjZUdRLDAuMjc3Nzc3OGVtRicvJSdyc3BhY2VHRlUtRiM2Ji1GLDYlUSJmRidGOEY7LUY/Ni1RMCZBcHBseUZ1bmN0aW9uO0YnRkJGREZHRklGS0ZNRk9GUS9GVFEmMC4wZW1GJy9GV0Zbby1JKG1mZW5jZWRHRiQ2JC1GIzYkLUYsNiVRInRGJ0Y4RjtGQkZCRkJGK0ZCLUY/Ni1RIixGJ0ZCRkQvRkhGOkZJRktGTUZPRlFGam4vRldRLDAuMzMzMzMzM2VtRictRiM2Jy1GLDYlUSJ5RidGOEY7Rj4tRiM2Ji1GLDYlUSJnRidGOEY7RmduRl1vRkJGK0ZCRitGQkYrRkJGK0ZC. For this you can use plot([f(t), g(t), t=a..b]); For example: plot([cos(t),sin(t),t=0..2*Pi]); In three dimensions an "implicit" curve would have to be given by two equations in the variables x,y,z, and intersectplot can be used for this: intersectplot(equation1, equation2, x=a1 .. b1, y=a2 .. b2, z = a3 .. b3); You can also have a parametric curve in 3 dimensions. For this you can use spacecurve in the plots package. spacecurve([x(t),y(t),z(t)], t=a..b); For example: spacecurve([sin(t),cos(t),sin(2*t)],t=0..2*Pi, colour=green,axes=box);

If you start Newton's method at a point close enough to a solution, it should converge rapidly to that solution. But how close? We investigated that for Newton's method in one variable. What about several variables? I'll look at a two-variable example, because the pictures will be easier to draw. F:= X -> <X[1]^2 * X[2] + X[1]*X[2]^2 - 2, X[1] + sin(X[2]) - 2>; J:= unapply(Jacobian(F(X),[X[1],X[2]]),X); W[0]:= <1, 2>; for count from 1 to 10 do W[count]:= Newt(W[count-1]) end do; Here are the curves corresponding to the two equations. with(plots): implicitplot([F(<x,y>)[1], F(<x,y>)[2]], x = 0 .. 4, y = -4 .. 4, colour = [red,blue]); Clearly there are two solutions in this rectangle, of which we've found one with Newt. I'll store that and look for the other one. Sol1:= W[10]; W[0]:= <2,-3>; for count from 1 to 10 do W[count]:= Newt(W[count-1]) end do; Sol2:=W[10]; Sol1,Sol2; Could there be more outside this rectangle? Take another look at the equations: F(<x,y>) = <0,0>; Since LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2JkYrLUYjNiYtRiw2JVEkc2luRicvJSdpdGFsaWNHUSZmYWxzZUYnLyUsbWF0aHZhcmlhbnRHUSdub3JtYWxGJy1JI21vR0YkNi1RMCZBcHBseUZ1bmN0aW9uO0YnRjkvJSZmZW5jZUdGOC8lKnNlcGFyYXRvckdGOC8lKXN0cmV0Y2h5R0Y4LyUqc3ltbWV0cmljR0Y4LyUobGFyZ2VvcEdGOC8lLm1vdmFibGVsaW1pdHNHRjgvJSdhY2NlbnRHRjgvJSdsc3BhY2VHUSYwLjBlbUYnLyUncnNwYWNlR0ZQLUkobWZlbmNlZEdGJDYkLUYjNiQtRiw2JVEieUYnL0Y3USV0cnVlRicvRjpRJ2l0YWxpY0YnRjlGOUY5RitGOUYrRjk= is always between -1 and 1, it's clear from the second equation that x must be between 1 and 3. Now the first equation is a quadratic in y, so any vertical line can only intersect its curve in at most two points. But for x between 1 and 3, we see both of those two points in our graph (one on each of the red curves). So there can't be any other solutions. Now to investigate our question: how close is "close enough"? One way to proceed is to plot circles of various radii centred at the solution together with their images under Newt. Suppose we find that for every radius LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEickYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy9GM1Enbm9ybWFsRic= with 0 < LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2Ji1GLDYlUSJyRicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUkjbW9HRiQ2LVElJmxlO0YnL0Y4USdub3JtYWxGJy8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGQi8lKXN0cmV0Y2h5R0ZCLyUqc3ltbWV0cmljR0ZCLyUobGFyZ2VvcEdGQi8lLm1vdmFibGVsaW1pdHNHRkIvJSdhY2NlbnRHRkIvJSdsc3BhY2VHUSwwLjI3Nzc3NzhlbUYnLyUncnNwYWNlR0ZRLUYsNiVRIlJGJ0Y0RjdGPkYrRj4=, the image of the circle is always inside a slightly smaller circle. Then, starting from any point whose distance from the solution is less than LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEiUkYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy9GM1Enbm9ybWFsRic=, you will always get convergence to the solution. I'll plot two parametric curves. One (in red) is the circle of radius R = 0.1 centred at the solution we found, the point Sol1, and the other (in blue) is the image of that under Newt. R:= 0.1: plot([[Sol1[1] + R*cos(t), Sol1[2] + R*sin(t), t = 0 .. 2*Pi], [ Newt(Sol1+<R*cos(t),R*sin(t)>)[1], Newt(Sol1+<R*cos(t),R*sin(t)>)[2], t = 0 .. 2*Pi]], colour=[red,blue]); The blue curve is the image of the red circle under Newt, and it's well within it.

LinearAlgebra package VectorCalculus package Jacobian (in VectorCalculus package) Determinant (in LinearAlgebra package) implicitplot3d (in plots package) intersectplot (in plots package) pointplot3d (in plots package) parametric version of plot spacecurve (in plots package)

LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkjbWlHRiQ2I1EhRic=