6元連立非線形方程式を解くプログラム

[firesick]

まず、2元連立非線形方程式
x^2+y^2-2=0
x^2-y^2-1=0
をx0=y0=1から始めて、連立Newton法で解くプログラムは、Jacobi行列の計算を数値微分で行い、

コード:

#include <iostream>
using namespace std;
#include <cstdio>
#include <cmath>

const int N=2;
const double EPS=1.0e-8;
const double h=1.0e-9;
const int MAX=40;

void Snewton(double *);
double fi(int, double *);
void Jacobi(double [][N],double *);
void LU(double *,double [][N],double *);

int main(){
	double x[N]={1,1};
	
	Snewton(x);
	
	printf("x=%f, y=%f \n", x[0],x[1]);
	return 0;
}

void Snewton(double x[]){
	int n,i,flag;
	double a[N][N],b[N],dx[N];
	printf("n      x          y\n0");
	printf(" %10.7f %10.7f\n",x[0],x[1]);
	
	for(n=1;n<MAX;n++){
	Jacobi(a,x);
	for(i=0;i<N;i++) b[i]=-fi(i,x);
	LU(dx,a,b);
	printf("%d",n);
	for(i=0;i<N;i++) x[i]+=dx[i];
	printf(" %10.7f %10.7f\n",x[0],x[1]);
	
	if(fabs(dx[0]) <= EPS && fabs(dx[1]) <= EPS) return;
}
}
	
double fi(int i,double x[]){
	switch(i){
	case 0: return x[0]*x[0]+x[1]*x[1]-2;
	case 1: return x[0]*x[0]-x[1]*x[1]-1;
	}
}

void Jacobi(double a[][N], double x[]){
	double xplush[N];	
	int i,j,k;
	
	for(i=0;i<N;i++){
	 for(j=0;j<N;j++){
	  for(k=0;k<N;k++) xplush[k]=x[k];
	  xplush[j]+=h;
	  a[i][j]=(fi(i,xplush)-fi(i,x))/h;
	  }
	}
} 

void LU(double dx[],double J[][N],double b[]){
	int i,j,k;
	double l[N][N],u[N][N],y[N];
	double sigma;
	
	for(j=0;j<N;j++) u[0][j]=J[0][j];
	for(j=1;j<N;j++) l[j][0]=J[j][0]/u[0][0];
	for(k=1;k<N;k++){
		for(j=k;j<N;j++)
			for(i=0;i<k;i++) J[j][k]-=l[j][i]*u[i][k];
		u[k][k]=J[k][k];
		for(j=k+1;j<N;j++){
			u[k][j]=J[k][j];
			for(i=0;i<k;i++) u[k][j]-=l[k][i]*u[i][j];
		}
	for(j=k+1;j<N;j++)l[j][k]=J[j][k]/u[k][k];
	}
	
	for(i=0;i<N;i++){
		y[i]=b[i];
		for(j=0;j<i;j++) y[i]-=l[i][j]*y[j];
	}
	
	for(k=1;k<=N;k++){
		i=N-k;
		sigma=y[i];
		for(j=i+1;j<N;j++)sigma-=u[i][j]*dx[j];
		dx[i]=sigma/u[i][i];
	}
}

これを改良して、6元連立非線形方程式を解けるようにしたいのですが、どうすればいいでしょうか。
出発近似解はx0=(0.5,2,1,-1,-0.5,-1)というふうに決めていいとします。
よろしくお願いします。

[みけCAT]

「連立Newton法」でググったらこういうのがありました。
3 非線型連立方程式の解(多元の場合)
これによると、今のプログラムのNを6にし、適切な行列を設定して反復計算をすればいいと思います。

[usao]

>どうすればいいでしょうか。

「現在困っていることは何でしょう？」というのを具体的に書かれたほうが良いと思います．

オフトピック

本題じゃないけど
> x^2+y^2-2=0
> x^2-y^2-1=0
みたいな　導関数の導出が簡単な式　を扱っているのであれば
安易に数値微分とか言わずに　導関数をそのまま実装したほうが良いような．