原理引用:https://blog.csdn.net/shiguangrenran1/article/details/8096206 Clear["Global`*"] ex31 = detaX*ex11 + detaY*ex12 + fun1 == 0 ex32 = detaX*ex21 + detaY*ex22 + fun2 == 0 ex50 = Solve[{ex31, ex32}, {detaX, detaY}] (* Solve[{fun1[x,y]==0,fun2[x,y]==0}]//Simplify *) Clear["Global`*"] fun1[x_, y_] := 5 x^5 - x^3 - 3 y^4 + y^2 - x + 665.0; fun2[x_, y_] := 4 x^7 + 7 x^5 - 3 y^5 - y^2 - x - 1046.0; Reduce[{fun1[x, y] == 0, fun2[x, y] == 0}, {x, y}, Reals] // N FindRoot[{fun1[x, y] == 0, fun2[x, y] == 0}, {x, 0}, {y, 10}] (* 牛顿迭代法 求解 *) ex11 = D[fun1[x0, y0], x0] ex12 = D[fun1[x0, y0], y0] ex21 = D[fun2[x0, y0], x0] ex22 = D[fun2[x0, y0], y0] ex31 = detaX*ex11 + detaY*ex12 + fun1[x0, y0] == 0 ex32 = detaX*ex21 + detaY*ex22 + fun2[x0, y0] == 0 ex50 = Solve[{ex31, ex32}, {detaX, detaY}] ; (* 解得 detaX detaY*) detaX = detaX /. ex50 detaY = detaY /. ex50 x0 = 10.0; y0 = 10.0; For[time = 1, time <= 20, time++, x0 = x0 + detaX; y0 = y0 + detaY; Print[time, " : ", Flatten[x0], ", ", Flatten[y0]];]