Solve a nonlinear system of equations with C/C++

Question

My system of equations is like this:

(x-a1)^2 + (y-b1)^2 = c1
(x-a2)^2 + (y-b2)^2 = c2

I know it is simple using Matlab:

solve((x-a1)^2 + (y-b1)^2 - c1, (x-a2)^2 + (y-b2)^2 - c2)

But how to solve this problem using C/C++? I know a math library lapack, but is for linear equations. Any suggestions?

Answer 1

You probably don't need to be using symbolic math in Matlab in the first place. If your parameters are given in terms of floating-point values and you're ultimately interested in obtaining a floating-point value as a result, numeric methods should be your first choice (there are reasons to use symbolic approaches, but I'm not going to discuss this here). For your example, you need a multi-dimensional nonlinear root-solver. In Matlab, fsolve (documentation) performs this function. You should probably start by ensuring that you can solve your problem in Matlab using this function – here's an example:

a = [1 2];
b = [3 4];
c = [5 6];
f = @(x)[(x(1)-a(1)).^2 + (x(2)-b(1)).^2-c(1);...
         (x(1)-a(2)).^2 + (x(2)-b(2)).^2-c(2)];
x0 = [1;1]; % Initial guess
opts  = optimset('Display','iter')
[x,fx] = fsolve(f,x0,opts)

which returns a single root with coordinates at

x =

   2.811249499748931
   1.688750500251424

Note that for these parameters, this system has multiple roots, but numeric methods generally only return one root. By adjusting the initial guess (here x0) or using using more advanced constrained optimization, the other roots can be found.

Matlab's fsolve uses a trust region method by default, the Trust-Region Dogleg method. In some cases the solver automatically switches algorithms with a warning message an you can specify the algorithm to use via optimset. Find the algorithm that works best for your problem.

As for C/C++, there are a variety of options. The GSL has support for multi-dimensional root-finding. You might look for code that implements the same or similar algorithm as you're using in fsolve. For example, it appears that MINPACK implements a Powell dogleg algorithm similar to fsolve's default method. There are many others, such as PNL and the quite interesting-looking O₂scl (documentation).