今回は x*x+y*y=2 と y=x の２つの方程式があった場合の解を、２変数のニュートン法で求める。

#include <stdio.h>
#include <math.h>

#define eps 1.0e-10             //精度
#define max 1000                //最大繰り返し回数

double f1(double x1,double x2)
{
        return x1*x1+x2*x2-2.0;
}

double f2(double x1,double x2)
{
        return x2-x1;
}

double J11(double x1,double x2)
{
        return 2.0*x1;
}

double J12(double x1,double x2)
{
        return 2.0*x2;
}

double J21(double x1,double x2)
{
        return -1.0;
}
double J22(double x1,double x2)
{
        return 1.0;
}

void newton2()
{
        int n;
        double x1,x2,c1,c2,det;

        n=0;
        x1=x2=5.0;

        do {
                if(n>max) {
                        printf("解が求まりませんでした。\n");
                        break; //無限ループを防ぐために一定回数繰り返したら終了
                }
                n++;

                det=J11(x1,x2)*J22(x1,x2)-J12(x1,x2)*J21(x1,x2);
                c1=-(J22(x1,x2)*f1(x1,x2)-J12(x1,x2)*f2(x1,x2))/det;
                c2=-(-J21(x1,x2)*f1(x1,x2)-J11(x1,x2)*f2(x1,x2))/det;
                x1+=c1;
                x2+=c2;

        } while (c1*c1+c2*c2 > eps );

        printf("解は %f と %f \n",x1,x2);
}

int main()
{
        newton2();

        return 0;
}

解は 1.000000 と 1.000000

ちなみに、x1 と x2 の初期値をうまく設定すれば、解 -1 と -1 も求まる。