-1

was doing a code in maple for non-linear systems using newtons I struggled a bit with this so my friend showed me his solution. I did exactly the same and I keep getting this error. I apologize for the long code, but I figured I'd add everything needed. If needed I can shorten it. I am just trying to do 5 iterations of this method, and then this error says "(in J)" and I am not sure what they mean by the 2nd argument

restart;

g := (x, y, z) -> 3*x - cos(y*z) - 1/2;

h := (x, y, z) -> x^2 - 81*(y + 0.1)^2 + sin(z) + 1.06;

i := (x, y, z) -> exp(-y*x) + 20*z + 10/3*Pi - 1;

A := (x, y, z) -> Vector[column](3, [g(x, y, z), h(x, y, z), i(x, y, z)]);


B := (x, y, z) -> Vector[column](3, [x, y, z]);

J := Matrix(3, 3, [[diff(g(x, y, z), x), diff(g(x, y, z), y), diff(g(x, y, z), z)], [diff(h(x, y, z), x), diff(h(x, y, z), y), diff(h(x, y, z), z)], [diff(i(x, y, z), x), diff(i(x, y, z), y), diff(i(x, y, z), z)]]);

J := (x, y, z) -> Matrix(3, 3, [[diff(g(x, y, z), x), diff(g(x, y, z), y), diff(g(x, y, z), z)], [diff(h(x, y, z), x), diff(h(x, y, z), y), diff(h(x, y, z), z)], [diff(i(x, y, z), x), diff(i(x, y, z), y), diff(i(x, y, z), z)]]);

C := (x, y, z) -> LinearAlgebra[MatrixInverse](J(x, y, z));

F := (x, y, z) -> evalf(B(x, y, z) - ((C(x, y, z)) . (A(x, y, z))));

x[0] := 0.5;
y[0] := 0.5;
z[0] := -0.5;

x[1] := F(x[0], y[0], z[0]);
Error, (in J) invalid input: diff received .5, which is not valid for its 2nd argument
acer
  • 6,671
  • 15
  • 15
Jos
  • 127
  • 6

1 Answers1

1

You could also use the Jacobian command from the VectorCalculus package, along with purely Vector form. But I'll try to make minimal changes to get your original to behave (eg. I'll leave all the unnecessary operator definitions).

The following works in my Maple 2021.0 if I Copy&Paste in as either 1D plaintext in a Worksheet or 2D Math in a Document.

restart;
g := (x,y,z) -> 3*x - cos(y*z) - 1/2:
h := (x,y,z) -> x^2 - 81*(y + 0.1)^2 + sin(z) + 1.06:
i := (x,y,z) -> exp(-y*x) + 20*z + 10/3*Pi - 1:

A := (x,y,z) -> Vector(3,[g(x,y,z),h(x,y,z),i(x,y,z)]):
B := (x,y,z) -> Vector(3,[x,y,z]):

J:=unapply(Matrix(3,3,
     [[diff(g(x,y,z),x),diff(g(x,y,z),y),diff(g(x,y,z),z)],
      [diff(h(x,y,z),x),diff(h(x,y,z),y),diff(h(x,y,z),z)],
      [diff(i(x,y,z),x),diff(i(x,y,z),y),diff(i(x,y,z),z)]]),
       [x,y,z]):

C := (x,y,z) -> LinearAlgebra:-MatrixInverse(J(x,y,z)):
F := (x,y,z) -> evalf(B(x,y,z) - ((C(x,y,z)) . (A(x,y,z)))):

x[0] := 0.5: y[0] := 0.5: z[0] := -0.5:
V[0] := Vector([x[0],y[0],z[0]]);

for ii from 1 to 7 do
  V[ii] := F(V[ii-1][1], V[ii-1][2], V[ii-1][3]);
  g(V[ii][1], V[ii][2], V[ii][3]),
  h(V[ii][1], V[ii][2], V[ii][3]),
  i(V[ii][1], V[ii][2], V[ii][3]);
end do;
acer
  • 6,671
  • 15
  • 15