Solve $f(x)=0$ for $x$.
Solve this problem using a sequence of $s$ values, starting at $s=0$.
Again, we have: $$\phi(x,s) = sf(x) + (1-s)g(x),$$
Take the derivative: $$d\phi = \frac{\partial\phi}{\partial x}dx + \frac{\partial\phi}{\partial s}ds = 0.$$
Divide by $ds$ and solve for $dx/ds$: $$\frac{dx}{ds} = -\left(\frac{\partial\phi}{\partial x}\right)^{-1}\left(\frac{\partial\phi}{\partial s}\right)_x.$$
If we consider $x$, $f(x)$ and $g(x)$ as vectors, then we have $$\frac{dx}{ds} = -J^{-1}\left(\frac{\partial\phi}{\partial s}\right)_x,$$ where $J$ is the Jacobian matrix with elements $J_{i,j} = \partial \phi_i/\partial x_j.$
This gives $$\frac{dx}{ds} = -J^{-1}\left(f(x)-g(x)\right).$$
We then solve the ODE (system) from $s=0$ to $s=1$ with initial condition $x_0$.
Use the following for $g(x)$, which has solution $x_0$, which is the initial condition for $x$ when $s=0$:
$$g(x) = f(x)-f(x_0).$$In this case $$\phi = f(x) + (s-1)f(x_0)$$
And we have
$$\frac{dx}{ds} = -J^{-1}f(x_0).$$with $J_{i,j} = \partial f_i/\partial x_j.$
Solve the following three equations in three unknowns using Newton, Sympy's fsolve, and Approch 2 homotopy-continuation method.
$$x^2 + y^2 = 1,$$$$xy + yz = -1.1,$$$$y^2 + z^2 = 2.$$using Plots
using LinearAlgebra
using NLsolve
using DifferentialEquations
#--------------- Analytic Jacobian
function jaca(f,xyz)
x,y,z = xyz
J = [2*x 2*y 0;
y x+z y;
0 2*y 2*z]
return J
end
#--------------- Numerical Jacobian
function jacn(f, x)
eps = 1E-8
n = length(x)
J = zeros(n,n)
for j in 1:n
h = eps * abs(x[j])
dx = zeros(n)
dx[j] = h
J[:,j] = (f(x.+dx) .- f(x)) / h
end
return J
end
#------------------ Simple Newton Solver
function newton(f, x0, tol=1E-6)
maxit = 50 # don't do more than this many iterations
x = x0 # initialize the solution x
for nit in 1:maxit
J = jaca(f, x)
dx = J\-f(x)
xnew = x .+ dx
err = norm((xnew.-x)./xnew)
x = xnew
if err < tol
break
end
if nit == maxit
println("warning, not converged in $maxit iterations")
end
end
return x
end
#------------------- Define the system to solve: f(x)
function f1f2f3(xyz)
x = xyz[1]
y = xyz[2]
z = xyz[3]
f1 = x*x + y*y - 1.
f2 = x*y + y*z + 1.1
f3 = y*y + z*z - 2.
return [f1, f2, f3]
end
#------------------ Right hand side function of the ODE to solve
function rhsf(x,t, f,f0)
J = jaca(f,x)
return -inv(J)*f0
return -J\f0
end
#------------------ Solve with Newton and fsolve
xyz0 = [1,2,3.] # initial guess
xyz0 = [-1., 2, -2]
xyz0 = [10., 2, 2]
xyz_nt = newton(f1f2f3, xyz0)
xyz_py = nlsolve( (eqn, xyz) -> eqn[:]=f1f2f3(xyz), xyz0).zero
#------------------ Homotopy solve
f0 = f1f2f3(xyz0)
s = [0,1.]
prob = ODEProblem((xyz, p, t)->rhsf(xyz,t,f1f2f3,f0), xyz0, s)
sol = solve(prob)
xyz_hc = sol.u[end]
println(xyz_nt)
println(xyz_py)
print(xyz_hc)
[0.8641393541884639, -0.5032525971544951, 1.3216417152380049] [0.8641393541892285, -0.5032525971540107, 1.3216417152383118] [0.8650474959675523, -0.5027938324351386, 1.32188193536268]
f1f2f3(xyz_hc)
3-element Vector{Float64}: 0.0011088082145467304 0.00042536995627417973 0.00017348897299873656