静岡理工科大学 菅沼ホーム 目次 索引

非線形方程式(ニュートン法)

    1. A. C++
    2. B. Java
    3. C. JavaScript
    4. D. PHP
    5. E. Ruby
    6. F. Python
    7. G. C#
    8. H. VB

  プログラムは,f(x) = exp(x) - 3x = 0 の根をニュートン法で求めた例です.多次元の場合に対するプログラムは,3 点 ( 0.5, 1.0 ),( 0.0, 1.5 ),( 0.5, 2.0 ) を通る円の中心座標と半径を多次元のニュートン法で求めた例です.

  1. C++

    /****************************/
    /* ニュートン法             */
    /* (exp(x)-3.0*x=0の根)   */
    /*      coded by Y.Suganuma */
    /****************************/
    #include <stdio.h>
    
    double newton(double(*)(double), double(*)(double), double, double,
                  double, int, int *);
    double snx(double);
    double dsnx(double);
    
    int main()
    {
    	double eps1, eps2, x, x0;
    	int max, ind;
    
    	eps1 = 1.0e-7;
    	eps2 = 1.0e-10;
    	max  = 20;
    	x0   = 0.0;
    
    	x = newton(snx, dsnx, x0, eps1, eps2, max, &ind);
    
    	printf("ind=%d  x=%f  f=%f  df=%f\n", ind, x, snx(x), dsnx(x));
    
    	return 0;
    }
    
    /*****************************************************/
    /* Newton法による非線形方程式(f(x)=0)の解            */
    /*      f : f(x)を計算する関数名                     */
    /*      df : f(x)の微分を計算する関数名              */
    /*      x0 : 初期値                                  */
    /*      eps1 : 終了条件1(|x(k+1)-x(k)|<eps1)   */
    /*      eps2 : 終了条件2(|f(x(k))|<eps2)       */
    /*      max : 最大試行回数                           */
    /*      ind : 実際の試行回数                         */
    /*            (負の時は解を得ることができなかった) */
    /*      return : 解                                  */
    /*****************************************************/
    #include <math.h>
    
    double newton(double(*f)(double), double(*df)(double), double x0,
                  double eps1, double eps2, int max, int *ind)
    {
    	double g, dg, x, x1;
    	int sw;
    
    	x1   = x0;
    	x    = x1;
    	*ind = 0;
    	sw   = 0;
    
    	while (sw == 0 && *ind >= 0) {
    
    		sw    = 1;
    		*ind += 1;
    		g     = (*f)(x1);
    
    		if (fabs(g) > eps2) {
    			if (*ind <= max) {
    				dg = (*df)(x1);
    				if (fabs(dg) > eps2) {
    					x = x1 - g / dg;
    					if (fabs(x-x1) > eps1 && fabs(x-x1) > eps1*fabs(x)) {
    						x1 = x;
    						sw = 0;
    					}
    				}
    				else
    					*ind = -1;
    			}
    			else
    				*ind = -1;
    		}
    	}
    
    	return x;
    }
    
    /************************/
    /* 関数値(f(x))の計算 */
    /************************/
    double snx(double x)
    {
    	double y;
    	y = exp(x) - 3.0 * x;
    	return y;
    }
    
    /********************/
    /* 関数の微分の計算 */
    /********************/
    double dsnx(double x)
    {
    	double y;
    	y = exp(x) - 3.0;
    	return y;
    }
    			
    多次元の場合
    /***************************************/
    /* 多次元ニュートン法                  */
    /* (3点(0.5,1.0),(0.0,1.5),(0.5,2.0) */
    /*   を通る円の中心と半径)            */
    /*      coded by Y.Suganuma            */
    /***************************************/
    #include <stdio.h>
    
    int newton(void(*)(double *, double *), int(*)(double *, double **, double),
               double *, double, double, int, double **, int, double *);
    void snx(double *, double *);
    int dsnx(double *, double **, double);
    int gauss(double **, int, int, double);
    
    int main()
    {
    	double eps1, eps2, x[3], x0[3], **w;
    	int i1, max, ind;
    
    	eps1  = 1.0e-7;
    	eps2  = 1.0e-10;
    	max   = 20;
    	x0[0] = 0.0;
    	x0[1] = 0.0;
    	x0[2] = 2.0;
    	w     = new double * [3];
    	for (i1 = 0; i1 < 3; i1++)
    		w[i1] = new double [6];
    
    	ind = newton(snx, dsnx, x0, eps1, eps2, max, w, 3, x);
    
    	printf("ind = %d\n", ind);
    	printf("x");
    	for (i1 = 0; i1 < 3; i1++)
    		printf(" %f", x[i1]);
    	printf("\nf ");
    	snx(x, x0);
    	for (i1 = 0; i1 < 3; i1++)
    		printf(" %f", x0[i1]);
    	printf("\n");
    
    	for (i1 = 0; i1 < 3; i1++)
    		delete [] w[i1];
    	delete [] w;
    
    	return 0;
    }
    
    /*****************************************************/
    /* Newton法による非線形方程式(f(x)=0)の解            */
    /*      f : f(x)を計算する関数名                     */
    /*      df : f(x)の微分を計算する関数名              */
    /*      x0 : 初期値                                  */
    /*      eps1 : 終了条件1(|x(k+1)-x(k)|<eps1)   */
    /*      eps2 : 終了条件2(|f(x(k))|<eps2)       */
    /*      max : 最大試行回数                           */
    /*      w : f(x)の微分の逆行列を入れる (n x 2n)      */
    /*      n : xの次数                                  */
    /*      x : 解                                       */
    /*      return : 実際の試行回数                      */
    /*            (負の時は解を得ることができなかった) */
    /*****************************************************/
    #include <math.h>
    
    int newton(void(*f)(double *, double *), int(*df)(double *, double **, double),
               double *x0, double eps1, double eps2, int max, double **w,
               int n, double *x)
    {
    	double *g, *x1;
    	int i1, i2, ind = 0, sw, sw1;
    
    	g  = new double [n];
    	x1 = new double [n];
    	sw = 0;
    	for (i1 = 0; i1 < n; i1++) {
    		x1[i1] = x0[i1];
    		x[i1]  = x1[i1];
    	}
    
    	while (sw == 0 && ind >= 0) {
    
    		sw = 1;
    		ind++;
    		(*f)(x1, g);
    
    		sw1 = 0;
    		for (i1 = 0; i1 < n && sw1 == 0; i1++) {
    			if (fabs(g[i1]) > eps2)
    				sw1 = 1;
    		}
    		if (sw1 > 0) {
    			if (ind <= max) {
    				sw1 = (*df)(x1, w, eps2);
    				if (sw1 == 0) {
    					for (i1 = 0; i1 < n; i1++) {
    						x[i1] = x1[i1];
    						for (i2 = 0; i2 < n; i2++)
    							x[i1] -= w[i1][i2+n] * g[i2];
    					}
    					sw1 = 0;
    					for (i1 = 0; i1 < n && sw1 == 0; i1++) {
    						if (fabs(x[i1]-x1[i1]) > eps1 && fabs(x[i1]-x1[i1]) > eps1*fabs(x[i1]))
    							sw1 = 1;
    					}
    					if (sw1 > 0) {
    						for (i1 = 0; i1 < n; i1++)
    							x1[i1] = x[i1];
    						sw = 0;
    					}
    				}
    				else
    					ind = -1;
    			}
    			else
    				ind = -1;
    		}
    	}
    
    	delete [] g;
    	delete [] x1;
    
    	return ind;
    }
    
    /************************/
    /* 関数値(f(x))の計算 */
    /************************/
    void snx(double *x, double *y)
    {
    	y[0] = (0.5 - x[0]) * (0.5 - x[0]) + (1.0 - x[1]) * (1.0 - x[1]) - x[2] * x[2];
    	y[1] = (0.0 - x[0]) * (0.0 - x[0]) + (1.5 - x[1]) * (1.5 - x[1]) - x[2] * x[2];
    	y[2] = (0.5 - x[0]) * (0.5 - x[0]) + (2.0 - x[1]) * (2.0 - x[1]) - x[2] * x[2];
    }
    
    /*****************************************/
    /* 関数の微分の計算                      */
    /*      x : 変数値                       */
    /*      w : 微分の逆行列(後半)         */
    /*      eps : 逆行列の存在を判定する規準 */
    /*      return : =0 : 正常               */
    /*               =1 : 逆行列が存在しない */
    /*****************************************/
    int dsnx(double *x, double **w, double eps)
    {
    	int i1, i2, sw;
    
    	for (i1 = 0; i1 < 3; i1++) {
    		for (i2 = 0; i2 < 3; i2++)
    			w[i1][i2+3] = 0.0;
    		w[i1][i1+3] = 1.0;
    	}
    
    	w[0][0] = -2.0 * (0.5 - x[0]);
    	w[0][1] = -2.0 * (1.0 - x[1]);
    	w[0][2] = -2.0 * x[2];
    	w[1][0] = -2.0 * (0.0 - x[0]);
    	w[1][1] = -2.0 * (1.5 - x[1]);
    	w[1][2] = -2.0 * x[2];
    	w[2][0] = -2.0 * (0.5 - x[0]);
    	w[2][1] = -2.0 * (2.0 - x[1]);
    	w[2][2] = -2.0 * x[2];
    
    	sw = gauss(w, 3, 3, eps);
    
    	return sw;
    }
    
    /*******************************************************/
    /* 線形連立方程式を解く(逆行列を求める)              */
    /*      w : 方程式の左辺及び右辺                       */
    /*      n : 方程式の数                                 */
    /*      m : 方程式の右辺の列の数                       */
    /*      eps : 逆行列の存在を判定する規準               */
    /*      return : =0 : 正常                             */
    /*               =1 : 逆行列が存在しない               */
    /*******************************************************/
    #include <math.h>
    
    int gauss(double **w, int n, int m, double eps)
    {
    	double y1, y2;
    	int ind = 0, nm, m1, m2, i1, i2, i3;
    
    	nm = n + m;
    
    	for (i1 = 0; i1 < n && ind == 0; i1++) {
    
    		y1 = .0;
    		m1 = i1 + 1;
    		m2 = 0;
    					// ピボット要素の選択
    		for (i2 = i1; i2 < n; i2++) {
    			y2 = fabs(w[i2][i1]);
    			if (y1 < y2) {
    				y1 = y2;
    				m2 = i2;
    			}
    		}
    					// 逆行列が存在しない
    		if (y1 < eps)
    			ind = 1;
    					// 逆行列が存在する
    		else {
    						// 行の入れ替え
    			for (i2 = i1; i2 < nm; i2++) {
    				y1        = w[i1][i2];
    				w[i1][i2] = w[m2][i2];
    				w[m2][i2] = y1;
    			}
    						// 掃き出し操作
    			y1 = 1.0 / w[i1][i1];
    
    			for (i2 = m1; i2 < nm; i2++)
    				w[i1][i2] *= y1;
    
    			for (i2 = 0; i2 < n; i2++) {
    				if (i2 != i1) {
    					for (i3 = m1; i3 < nm; i3++)
    						w[i2][i3] -= w[i2][i1] * w[i1][i3];
    				}
    			}
    		}
    	}
    
    	return(ind);
    }
    			

  2. Java

    /****************************/
    /* ニュートン法             */
    /* (exp(x)-3.0*x=0の根)   */
    /*      coded by Y.Suganuma */
    /****************************/
    import java.io.*;
    
    public class Test {
    	public static void main(String args[]) throws IOException
    	{
    		double eps1, eps2, x, x0;
    		int max, ind[] = new int [1];
    
    		eps1 = 1.0e-7;
    		eps2 = 1.0e-10;
    		max  = 30;
    		x0   = 0.0;
    					// 関数値を計算するクラス
    		Kansu kn = new Kansu();
    					// ニュートン法の実行
    		x = kn.newton(x0, eps1, eps2, max, ind);
    					// 出力
    		System.out.println("ind=" + ind[0] + "  x=" + x + "  f=" + kn.snx(x) + "  df=" + kn.dsnx(x));
    	}
    }
    
    /******************************/
    /* 関数値およびその微分の計算 */
    /******************************/
    class Kansu extends Newton
    {
    					// 関数値(f(x))の計算
    	double snx(double x)
    	{
    		double y = Math.exp(x) - 3.0 * x;
    		return y;
    	}
    					// 関数の微分の計算
    	double dsnx(double x)
    	{
    		double y = Math.exp(x) - 3.0;
    		return y;
    	}
    }
    
    abstract class Newton
    {
    	/*****************************************************/
    	/* Newton法による非線形方程式(f(x)=0)の解            */
    	/*      x1 : 初期値                                  */
    	/*      eps1 : 終了条件1(|x(k+1)-x(k)|<eps1)   */
    	/*      eps2 : 終了条件2(|f(x(k))|<eps2)       */
    	/*      max : 最大試行回数                           */
    	/*      ind : 実際の試行回数                         */
    	/*            (負の時は解を得ることができなかった) */
    	/*      return : 解                                  */
    	/*****************************************************/
    
    	abstract double snx(double x);
    	abstract double dsnx(double x);
    
    	double newton(double x1, double eps1, double eps2, int max, int ind[])
    	{
    		double g, dg, x;
    		int sw;
    
    		x      = x1;
    		ind[0] = 0;
    		sw     = 0;
    
    		while (sw == 0 && ind[0] >= 0) {
    
    			ind[0]++;
    			sw = 1;
    			g  = snx(x1);
    
    			if (Math.abs(g) > eps2) {
    				if (ind[0] <= max) {
    					dg = dsnx(x1);
    					if (Math.abs(dg) > eps2) {
    						x = x1 - g / dg;
    						if (Math.abs(x-x1) > eps1 && Math.abs(x-x1) > eps1*Math.abs(x)) {
    							x1 = x;
    							sw = 0;
    						}
    					}
    					else
    						ind[0] = -1;
    				}
    				else
    					ind[0] = -1;
    			}
    		}
    
    		return x;
    	}
    }
    			
    多次元の場合
    /***************************************/
    /* 多次元ニュートン法                  */
    /* (3点(0.5,1.0),(0.0,1.5),(0.5,2.0) */
    /*   を通る円の中心と半径)            */
    /*      coded by Y.Suganuma            */
    /***************************************/
    import java.io.*;
    
    public class Test_m {
    	public static void main(String args[]) throws IOException
    	{
    		double eps1, eps2, x[], x0[], w[][];
    		int i1, max, ind;
    
    		x     = new double [3];
    		x0    = new double [3];
    		w     = new double [3][6];
    		eps1  = 1.0e-7;
    		eps2  = 1.0e-10;
    		max   = 20;
    		x0[0] = 0.0;
    		x0[1] = 0.0;
    		x0[2] = 2.0;
    
    		Kansu kn = new Kansu();
    
    		ind = kn.newton_m(x0, eps1, eps2, max, w, 3, x);
    
    		System.out.println("ind = " + ind);
    		System.out.print("x");
    		for (i1 = 0; i1 < 3; i1++)
    			System.out.print(" " + x[i1]);
    		System.out.println("\n");
    		kn.snx(x, x0);
    		for (i1 = 0; i1 < 3; i1++)
    			System.out.print(" " + x0[i1]);
    		System.out.println();
    	}
    }
    
    /******************************/
    /* 関数値およびその微分の計算 */
    /******************************/
    class Kansu extends Newton
    {
    					// 関数値(f(x))の計算
    	void snx(double x[], double y[])
    	{
    		y[0] = (0.5 - x[0]) * (0.5 - x[0]) + (1.0 - x[1]) * (1.0 - x[1]) - x[2] * x[2];
    		y[1] = (0.0 - x[0]) * (0.0 - x[0]) + (1.5 - x[1]) * (1.5 - x[1]) - x[2] * x[2];
    		y[2] = (0.5 - x[0]) * (0.5 - x[0]) + (2.0 - x[1]) * (2.0 - x[1]) - x[2] * x[2];
    	}
    					// 関数の微分の計算
    	int dsnx(double x[], double w[][], double eps)
    	{
    		for (int i1 = 0; i1 < 3; i1++) {
    			for (int i2 = 0; i2 < 3; i2++)
    				w[i1][i2+3] = 0.0;
    			w[i1][i1+3] = 1.0;
    		}
    
    		w[0][0] = -2.0 * (0.5 - x[0]);
    		w[0][1] = -2.0 * (1.0 - x[1]);
    		w[0][2] = -2.0 * x[2];
    		w[1][0] = -2.0 * (0.0 - x[0]);
    		w[1][1] = -2.0 * (1.5 - x[1]);
    		w[1][2] = -2.0 * x[2];
    		w[2][0] = -2.0 * (0.5 - x[0]);
    		w[2][1] = -2.0 * (2.0 - x[1]);
    		w[2][2] = -2.0 * x[2];
    
    		int ind = gauss(w, 3, 3, eps);
    
    		return ind;
    	}
    }
    
    abstract class Newton
    {
    	/*****************************************************/
    	/* Newton法による多次元非線形方程式(f(x)=0)の解      */
    	/*      x1 : 初期値                                  */
    	/*      eps1 : 終了条件1(|x(k+1)-x(k)|<eps1)   */
    	/*      eps2 : 終了条件2(|f(x(k))|<eps2)       */
    	/*      max : 最大試行回数                           */
    	/*      kn1 : 関数を計算するクラスオブジェクト       */
    	/*      kn2 : 関数の微分を計算するクラスオブジェクト */
    	/*      w : f(x)の微分の逆行列を入れる (n x 2n)      */
    	/*      n : xの次数                                  */
    	/*      x : 解                                       */
    	/*      return : 実際の試行回数                      */
    	/*            (負の時は解を得ることができなかった) */
    	/*****************************************************/
    
    	abstract void snx(double x[], double y[]);
    	abstract int dsnx(double x[], double w[][], double eps);
    
    	int newton_m(double x0[], double eps1, double eps2, int max, double w[][], int n, double x[])
    	{
    		double g[], x1[];
    		int i1, i2, ind = 0, sw, sw1;
    
    		g  = new double [n];
    		x1 = new double [n];
    		sw = 0;
    		for (i1 = 0; i1 < n; i1++) {
    			x1[i1] = x0[i1];
    			x[i1]  = x1[i1];
    		}
    
    		while (sw == 0 && ind >= 0) {
    
    			sw = 1;
    			ind++;
    			snx(x1, g);
    
    			sw1 = 0;
    			for (i1 = 0; i1 < n && sw1 == 0; i1++) {
    				if (Math.abs(g[i1]) > eps2)
    					sw1 = 1;
    			}
    			if (sw1 > 0) {
    				if (ind <= max) {
    					sw1 = dsnx(x1, w, eps2);
    					if (sw1 == 0) {
    						for (i1 = 0; i1 < n; i1++) {
    							x[i1] = x1[i1];
    							for (i2 = 0; i2 < n; i2++)
    								x[i1] -= w[i1][i2+n] * g[i2];
    						}
    						sw1 = 0;
    						for (i1 = 0; i1 < n && sw1 == 0; i1++) {
    							if (Math.abs(x[i1]-x1[i1]) > eps1 &&
                                    Math.abs(x[i1]-x1[i1]) > eps1*Math.abs(x[i1]))
    								sw1 = 1;
    						}
    						if (sw1 > 0) {
    							for (i1 = 0; i1 < n; i1++)
    								x1[i1] = x[i1];
    							sw = 0;
    						}
    					}
    					else
    						ind = -1;
    				}
    				else
    					ind = -1;
    			}
    		}
    
    		return ind;
    	}
    
    	/*******************************************************/
    	/* 線形連立方程式を解く(逆行列を求める)              */
    	/*      w : 方程式の左辺及び右辺                       */
    	/*      n : 方程式の数                                 */
    	/*      m : 方程式の右辺の列の数                       */
    	/*      eps : 逆行列の存在を判定する規準               */
    	/*      return : =0 : 正常                             */
    	/*               =1 : 逆行列が存在しない               */
    	/*******************************************************/
    	int gauss(double w[][], int n, int m, double eps)
    	{
    		double y1, y2;
    		int ind = 0, nm, m1, m2, i1, i2, i3;
    
    		nm = n + m;
    
    		for (i1 = 0; i1 < n && ind == 0; i1++) {
    
    			y1 = .0;
    			m1 = i1 + 1;
    			m2 = 0;
    						// ピボット要素の選択
    			for (i2 = i1; i2 < n; i2++) {
    				y2 = Math.abs(w[i2][i1]);
    				if (y1 < y2) {
    					y1 = y2;
    					m2 = i2;
    				}
    			}
    						// 逆行列が存在しない
    			if (y1 < eps)
    				ind = 1;
    						// 逆行列が存在する
    			else {
    							// 行の入れ替え
    				for (i2 = i1; i2 < nm; i2++) {
    					y1        = w[i1][i2];
    					w[i1][i2] = w[m2][i2];
    					w[m2][i2] = y1;
    				}
    							// 掃き出し操作
    				y1 = 1.0 / w[i1][i1];
    
    				for (i2 = m1; i2 < nm; i2++)
    					w[i1][i2] *= y1;
    
    				for (i2 = 0; i2 < n; i2++) {
    					if (i2 != i1) {
    						for (i3 = m1; i3 < nm; i3++)
    							w[i2][i3] -= w[i2][i1] * w[i1][i3];
    					}
    				}
    			}
    		}
    
    		return(ind);
    	}
    }
    			

  3. JavaScript

      ここをクリック(多次元の場合は,ここをクリック)すると,JavaScript の仕様に適合した形で解を求めたい式を入力することによって,任意の非線形方程式の解を画面上で求めることができます.
    <!DOCTYPE HTML>
    
    <HTML>
    
    <HEAD>
    
    	<TITLE>非線形方程式(ニュートン法)</TITLE>
    	<META HTTP-EQUIV="Content-Type" CONTENT="text/html; charset=utf-8">
    	<SCRIPT TYPE="text/javascript">
    		str1 = "";
    		str2 = "";
    		ind = 0;   // 実際の試行回数(負の時は解を得ることができなかった)
    		function main()
    		{
    			str1 = document.getElementById("func").value + ";";
    			str2 = document.getElementById("dfunc").value + ";";
    					// データの設定
    			let eps1 = 1.0e-7;
    			let eps2 = 1.0e-10;
    			let max  = parseInt(document.getElementById("max").value);
    			let x0   = parseFloat(document.getElementById("x0").value);
    			ind      = 0;
    					// 実行
    			let x = newton(x0, eps1, eps2, max, snx, dsnx);
    					// 結果
    			if (ind < 0)
    				document.getElementById("res").value = "解を求めることができませんでした!";
    			else {
    				let c = 100000;
    				let s1 = Math.round(x * c) / c;
    				document.getElementById("res").value = " x = " + s1 + ", 収束回数:" + ind;
    			}
    		}
    
    		/****************/
    		/* 関数値の計算 */
    		/****************/
    		function snx(x)
    		{
    			let y = eval(str1);
    			return y;
    		}
    
    		/********************/
    		/* 関数の微分の計算 */
    		/********************/
    		function dsnx(x)
    		{
    			let y = eval(str2);
    			return y;
    		}
    
    		/***************************************************/
    		/* Newton法による非線形方程式(f(x)=0)の解          */
    		/*      x0 : 初期値                                */
    		/*      eps1 : 終了条件1(|x(k+1)-x(k)|<eps1) */
    		/*      eps2 : 終了条件2(|f(x(k))|<eps2)     */
    		/*      max : 最大試行回数                         */
    		/*      fn : 関数値を計算する関数         */
    		/*      dfn : 関数の微分値を計算する関数           */
    		/*      return : 解                                */
    		/***************************************************/
    		function newton(x0, eps1, eps2, max, fn, dfn)
    		{
    			let x1   = x0;
    			let x    = x1;
    			let sw   = 0;
    
    			while (sw == 0 && ind >= 0) {
    
    				sw   = 1;
    				ind += 1;
    				let g     = fn(x1);
    
    				if (Math.abs(g) > eps2) {
    					if (ind <= max) {
    						let dg = dfn(x1);
    						if (Math.abs(dg) > eps2) {
    							x = x1 - g / dg;
    							if (Math.abs(x-x1) > eps1 && Math.abs(x-x1) > eps1*Math.abs(x)) {
    								x1 = x;
    								sw = 0;
    							}
    						}
    						else
    							ind = -1;
    					}
    					else
    						ind = -1;
    				}
    			}
    
    			return x;
    		}
    	</SCRIPT>
    
    </HEAD>
    
    <BODY STYLE="font-size: 130%; background-color: #eeffee;">
    
    	<H2 STYLE="text-align:center"><B>非線形方程式(ニュートン法)</B></H2>
    
    	<DL>
    		<DT>  テキストフィールドには,例として,以下に示すような非線形方程式の解を求める場合に対する値が設定されています.他の問題を実行する場合は,それらを適切に修正してください.なお,式や式の微分は,JavaScript の仕様に適合した形式で記述してあることに注意してください.
    		<P STYLE="text-align:center"><IMG SRC="newton.gif"></P>
    	</DL>
    
    	<DIV STYLE="text-align:center">
    		初期値:<INPUT ID="x0" STYLE="font-size: 100%" TYPE="text" SIZE="5" VALUE="0"> 
    		最大繰り返し回数:<INPUT ID="max" STYLE="font-size: 100%" TYPE="text" SIZE="5" VALUE="100"><BR><BR>
    		式: f(x) = <INPUT ID="func" STYLE="font-size: 100%" TYPE="text" SIZE="50" VALUE="Math.exp(x) - 3.0 * x">= 0<BR>
    		式の微分: f'(x) = <INPUT ID="dfunc" STYLE="font-size: 100%" TYPE="text" SIZE="50" VALUE="Math.exp(x) - 3.0">  
    		<BUTTON STYLE="font-size: 100%; background-color: pink" onClick="main()">実行</BUTTON><BR><BR>
    		結果:<INPUT ID="res" STYLE="font-size: 100%" TYPE="text" SIZE="30">
    	</DIV>
    
    </BODY>
    
    </HTML>
    			
    多次元の場合
    <!DOCTYPE HTML>
    
    <HTML>
    
    <HEAD>
    
    	<TITLE>非線形方程式(ニュートン法,多次元)</TITLE>
    	<META HTTP-EQUIV="Content-Type" CONTENT="text/html; charset=utf-8">
    	<SCRIPT TYPE="text/javascript">
    		str1 = new Array();
    		str2 = new Array();
    		function main()
    		{
    					// データの設定
    			let eps1 = 1.0e-7;
    			let eps2 = 1.0e-10;
    			let n    = parseInt(document.getElementById("order").value);
    			let max  = parseInt(document.getElementById("max").value);
    			let x    = new Array();
    			let x0   = new Array();
    			x0 = document.getElementById("x0").value.split(/ {1,}/)
    			for (let i1 = 0; i1 < n; i1++)
    				x0[i1] = parseFloat(x0[i1]);
    			str1 = document.getElementById("func").value.split(/\n{1,}/)
    			for (let i1 = 0; i1 < n; i1++)
    				str1[i1] = str1[i1] + ";";
    			let str3 = document.getElementById("dfunc").value.split(/\n{1,}/)
    			let w = new Array();
    			for (let i1 = 0; i1 < n; i1++) {
    				w[i1] = new Array();
    				str2[i1] = new Array();
    				for (let i2 = 0; i2 < n; i2++)
    					str2[i1][i2] = str3[i1*n+i2] + ";";
    			}
    					// 実行
    			let ind = newton(x0, eps1, eps2, max, w, n, x, snx, dsnx);
    					// 結果
    			if (ind < 0)
    				document.getElementById("res").value = "解を求めることができませんでした!";
    			else {
    				let c = 100000;
    				let str = "収束回数:" + ind + "\n";
    				str += "  x =";
    				for (let i1 = 0; i1 < n; i1++)
    					str = str + " " + Math.round(x[i1] * c) / c;
    				snx(n, x, x0);
    				str += "\n  f(x) =";
    				for (let i1 = 0; i1 < n; i1++)
    					str = str + " " + Math.round(x0[i1] * c) / c;
    				document.getElementById("res").value = str;
    			}
    		}
    
    		/****************/
    		/* 関数値の計算 */
    		/****************/
    		function snx(n, x, y)
    		{
    			for (let i1 = 0; i1 < n; i1++)
    				y[i1] = eval(str1[i1]);
    		}
    
    		/********************/
    		/* 関数の微分の計算 */
    		/********************/
    		function dsnx(n, x, w, eps)
    		{
    			for (let i1 = 0; i1 < n; i1++) {
    				for (let i2 = 0; i2 < n; i2++)
    					w[i1][i2+n] = 0.0;
    				w[i1][i1+n] = 1.0;
    			}
    
    			for (let i1 = 0; i1 < n; i1++) {
    				for (let i2 = 0; i2 < n; i2++)
    					w[i1][i2] = eval(str2[i1][i2]);
    			}
    
    			let ind = gauss(w, n, n, eps);
    
    			return ind;
    		}
    
    		/*****************************************************/
    		/* Newton法による非線形方程式(f(x)=0)の解            */
    		/*      x1 : 初期値                                  */
    		/*      eps1 : 終了条件1(|x(k+1)-x(k)|<eps1)   */
    		/*      eps2 : 終了条件2(|f(x(k))|<eps2)       */
    		/*      max : 最大試行回数                           */
    		/*      w : f(x)の微分の逆行列を入れる (n x 2n)      */
    		/*      n : xの次数                                  */
    		/*      x : 解                                       */
    		/*      fn : 関数値を計算する関数           */
    		/*      dfn : 関数の微分値を計算する関数             */
    		/*      return : 実際の試行回数                      */
    		/*            (負の時は解を得ることができなかった) */
    		/*****************************************************/
    		function newton(x0, eps1, eps2, max, w, n, x, fn, dfn)
    		{
    			let g  = new Array();
    			let x1 = new Array();
    			let sw = 0;
    			let ind = 0;
    			for (let i1 = 0; i1 < n; i1++) {
    				x1[i1] = x0[i1];
    				x[i1]  = x1[i1];
    			}
    
    			while (sw == 0 && ind >= 0) {
    
    				sw = 1;
    				ind++;
    				fn(n, x1, g);
    
    				let sw1 = 0;
    				for (let i1 = 0; i1 < n && sw1 == 0; i1++) {
    					if (Math.abs(g[i1]) > eps2)
    						sw1 = 1;
    				}
    				if (sw1 > 0) {
    					if (ind <= max) {
    						sw1 = dfn(n, x1, w, eps2);
    						if (sw1 == 0) {
    							for (i1 = 0; i1 < n; i1++) {
    								x[i1] = x1[i1];
    								for (i2 = 0; i2 < n; i2++)
    									x[i1] -= w[i1][i2+n] * g[i2];
    							}
    							sw1 = 0;
    							for (i1 = 0; i1 < n && sw1 == 0; i1++) {
    								if (Math.abs(x[i1]-x1[i1]) > eps1 &&
                    	                Math.abs(x[i1]-x1[i1]) > eps1*Math.abs(x[i1]))
    									sw1 = 1;
    							}
    							if (sw1 > 0) {
    								for (i1 = 0; i1 < n; i1++)
    									x1[i1] = x[i1];
    								sw = 0;
    							}
    						}
    						else
    							ind = -1;
    					}
    					else
    						ind = -1;
    				}
    			}
    
    			return ind;
    		}
    
    		/******************************************/
    		/* 線形連立方程式を解く(逆行列を求める) */
    		/*      w : 方程式の左辺及び右辺          */
    		/*      n : 方程式の数                    */
    		/*      m : 方程式の右辺の列の数          */
    		/*      eps : 逆行列の存在を判定する規準  */
    		/*      return : =0 : 正常                */
    		/*               =1 : 逆行列が存在しない  */
    		/******************************************/
    		function gauss(w, n, m, eps)
    		{
    			let ind = 0;
    			let nm = n + m;
    
    			for (let i1 = 0; i1 < n && ind == 0; i1++) {
    
    				let y1 = .0;
    				let m1 = i1 + 1;
    				let m2 = 0;
    						// ピボット要素の選択
    				for (let i2 = i1; i2 < n; i2++) {
    					let y2 = Math.abs(w[i2][i1]);
    					if (y1 < y2) {
    						y1 = y2;
    						m2 = i2;
    					}
    				}
    						// 逆行列が存在しない
    				if (y1 < eps)
    					ind = 1;
    						// 逆行列が存在する
    				else {
    							// 行の入れ替え
    					for (let i2 = i1; i2 < nm; i2++) {
    						y1        = w[i1][i2];
    						w[i1][i2] = w[m2][i2];
    						w[m2][i2] = y1;
    					}
    							// 掃き出し操作
    					y1 = 1.0 / w[i1][i1];
    
    					for (let i2 = m1; i2 < nm; i2++)
    						w[i1][i2] *= y1;
    
    					for (let i2 = 0; i2 < n; i2++) {
    						if (i2 != i1) {
    							for (let i3 = m1; i3 < nm; i3++)
    								w[i2][i3] -= w[i2][i1] * w[i1][i3];
    						}
    					}
    				}
    			}
    
    			return ind;
    		}
    	</SCRIPT>
    
    </HEAD>
    
    <BODY STYLE="font-size: 130%; background-color: #eeffee;">
    
    	<H2 STYLE="text-align:center"><B>非線形方程式(ニュートン法,多次元)</B></H2>
    
    	<DL>
    		<DT>  テキストフィールド及びテキストエリアには,例として,3 点 ( 0.5, 1.0 ),( 0.0, 1.5 ),( 0.5, 2.0 ) を通る円の中心座標と半径を求める場合に対する値が設定されています(他の問題を実行する場合は,それらを適切に修正してください).円の中心座標を (x, y ),半径を r とすると,以下に示す 3 つ式が成立しますので,これらの式を同時に満たす x, y, r を求めれば良いことになります.なお,プログラム上では,これらの変数が,それぞれ,配列の x[0], x[1], x[2] に対応しています.
    		<P STYLE="text-align:center"><IMG SRC="newton_m.gif"></P>
    		<DT>  関数の微分は行列となり,変数 x, y, r を,それぞれ,x<SUB>1</SUB>, x<SUB>2</SUB>, x<SUB>3</SUB> としたとき,その微分は以下のようになります.
    		<P STYLE="text-align:center"><IMG SRC="newton_m_bibun.gif"></P>
    		<DT>これらの式を,1 行 1 列目,1 行 2 列目,・・・,3 行 3 列目の順で入力して下さい.また,この例に示したように,関数及びその微分における変数は配列の形で記述して下さい.なお,式や式の微分は,JavaScript の仕様に適合した形式で記述してあることに注意してください.
    	</DL>
    
    	<DIV STYLE="text-align:center">
    		次数:<INPUT ID="order" STYLE="font-size: 100%" TYPE="text" SIZE="2" VALUE="3"> 
    		初期値:<INPUT ID="x0" STYLE="font-size: 100%" TYPE="text" SIZE="15" VALUE="0.0 0.0 2.0"> 
    		最大繰り返し回数:<INPUT ID="max" STYLE="font-size: 100%" TYPE="text" SIZE="3" VALUE="100"><BR><BR>
    		式: f(x) = <TEXTAREA ID="func" STYLE="font-size: 110%" COLS="70" ROWS="10">
    (0.5 - x[0]) * (0.5 - x[0]) + (1.0 - x[1]) * (1.0 - x[1]) - x[2] * x[2]
    (0.0 - x[0]) * (0.0 - x[0]) + (1.5 - x[1]) * (1.5 - x[1]) - x[2] * x[2]
    (0.5 - x[0]) * (0.5 - x[0]) + (2.0 - x[1]) * (2.0 - x[1]) - x[2] * x[2]
    		</TEXTAREA> = 0<BR><BR>
    		式の微分: f'(x) = <TEXTAREA ID="dfunc" STYLE="font-size: 110%" COLS="70" ROWS="10">
    -2.0 * (0.5 - x[0])
    -2.0 * (1.0 - x[1])
    -2.0 * x[2]
    -2.0 * (0.0 - x[0])
    -2.0 * (1.5 - x[1])
    -2.0 * x[2]
    -2.0 * (0.5 - x[0])
    -2.0 * (2.0 - x[1])
    -2.0 * x[2]
    		</TEXTAREA><BR><BR>
    		<BUTTON STYLE="font-size: 100%; background-color: pink" onClick="main()">実行</BUTTON><BR><BR>
    		結果:<TEXTAREA ID="res" STYLE="font-size: 110%" COLS="70" ROWS="10"></TEXTAREA>
    	</DIV>
    
    </BODY>
    
    </HTML>
    			

  4. PHP

    <?php
    
    /****************************/
    /* ニュートン法             */
    /* (exp(x)-3.0*x=0の根)   */
    /*      coded by Y.Suganuma */
    /****************************/
    
    	$eps1 = 1.0e-7;
    	$eps2 = 1.0e-10;
    	$max  = 20;
    	$x0   = 0.0;
    
    	$x = newton("snx", "dsnx", $x0, $eps1, $eps2, $max, $ind);
    
    	printf("ind=%d  x=%f  f=%f  df=%f\n", $ind, $x, snx($x), dsnx($x));
    
    /************************/
    /* 関数値(f(x))の計算 */
    /************************/
    function snx($x)
    {
    	return exp($x) - 3.0 * $x;
    }
    
    /********************/
    /* 関数の微分の計算 */
    /********************/
    function dsnx($x)
    {
    	return exp($x) - 3.0;
    }
    
    /*****************************************************/
    /* Newton法による非線形方程式(f(x)=0)の解            */
    /*      f : f(x)を計算する関数名                     */
    /*      df : f(x)の微分を計算する関数名              */
    /*      x0 : 初期値                                  */
    /*      eps1 : 終了条件1(|x(k+1)-x(k)|<eps1)   */
    /*      eps2 : 終了条件2(|f(x(k))|<eps2)       */
    /*      max : 最大試行回数                           */
    /*      ind : 実際の試行回数                         */
    /*            (負の時は解を得ることができなかった) */
    /*      return : 解                                  */
    /*****************************************************/
    function newton($f, $df, $x0, $eps1, $eps2, $max, &$ind)
    {
    	$x1  = $x0;
    	$x   = $x1;
    	$ind = 0;
    	$sw  = 0;
    
    	while ($sw == 0 && $ind >= 0) {
    
    		$sw   = 1;
    		$ind += 1;
    		$g    = $f($x1);
    
    		if (abs($g) > $eps2) {
    			if ($ind <= $max) {
    				$dg = $df($x1);
    				if (abs($dg) > $eps2) {
    					$x = $x1 - $g / $dg;
    					if (abs($x-$x1) > $eps1 && abs($x-$x1) > $eps1*abs($x)) {
    						$x1 = $x;
    						$sw = 0;
    					}
    				}
    				else
    					$ind = -1;
    			}
    			else
    				$ind = -1;
    		}
    	}
    
    	return $x;
    }
    
    ?>
    			
    多次元の場合
    <?php
    
    /***************************************/
    /* 多次元ニュートン法                  */
    /* (3点(0.5,1.0),(0.0,1.5),(0.5,2.0) */
    /*   を通る円の中心と半径)            */
    /*      coded by Y.Suganuma            */
    /***************************************/
    
    	$eps1  = 1.0e-7;
    	$eps2  = 1.0e-10;
    	$max   = 20;
    	$x     = array();
    	$x0    = array();
    	$x0[0] = 0.0;
    	$x0[1] = 0.0;
    	$x0[2] = 2.0;
    	$w     = array();
    	for ($i1 = 0; $i1 < 3; $i1++)
    		$w[$i1] = array();
    
    	$ind = newton("snx", "dsnx", $x0, $eps1, $eps2, $max, $w, 3, $x);
    
    	printf("ind = %d\n", $ind);
    	printf("x");
    	for ($i1 = 0; $i1 < 3; $i1++)
    		printf(" %f", $x[$i1]);
    	printf("\nf ");
    	snx($x, $x0);
    	for ($i1 = 0; $i1 < 3; $i1++)
    		printf(" %f", $x0[$i1]);
    	printf("\n");
    
    /*****************************************************/
    /* Newton法による非線形方程式(f(x)=0)の解            */
    /*      f : f(x)を計算する関数名                     */
    /*      df : f(x)の微分を計算する関数名              */
    /*      x0 : 初期値                                  */
    /*      eps1 : 終了条件1(|x(k+1)-x(k)|<eps1)   */
    /*      eps2 : 終了条件2(|f(x(k))|<eps2)       */
    /*      max : 最大試行回数                           */
    /*      w : f(x)の微分の逆行列を入れる (n x 2n)      */
    /*      n : xの次数                                  */
    /*      x : 解                                       */
    /*      return : 実際の試行回数                      */
    /*            (負の時は解を得ることができなかった) */
    /*****************************************************/
    function newton($f, $df, $x0, $eps1, $eps2, $max, $w, $n, &$x)
    {
    	$ind = 0;
    	$g   = array($n);
    	$x1  = array($n);
    	$sw  = 0;
    	for ($i1 = 0; $i1 < $n; $i1++) {
    		$x1[$i1] = $x0[$i1];
    		$x[$i1]  = $x1[$i1];
    	}
    
    	while ($sw == 0 && $ind >= 0) {
    
    		$sw = 1;
    		$ind++;
    		$f($x1, $g);
    
    		$sw1 = 0;
    		for ($i1 = 0; $i1 < $n && $sw1 == 0; $i1++) {
    			if (abs($g[$i1]) > $eps2)
    				$sw1 = 1;
    		}
    		if ($sw1 > 0) {
    			if ($ind <= $max) {
    				$sw1 = $df($x1, $w, $eps2);
    				if ($sw1 == 0) {
    					for ($i1 = 0; $i1 < $n; $i1++) {
    						$x[$i1] = $x1[$i1];
    						for ($i2 = 0; $i2 < $n; $i2++)
    							$x[$i1] -= $w[$i1][$i2+$n] * $g[$i2];
    					}
    					$sw1 = 0;
    					for ($i1 = 0; $i1 < $n && $sw1 == 0; $i1++) {
    						if (abs($x[$i1]-$x1[$i1]) > $eps1 && abs($x[$i1]-$x1[$i1]) > $eps1*abs($x[$i1]))
    							$sw1 = 1;
    					}
    					if ($sw1 > 0) {
    						for ($i1 = 0; $i1 < $n; $i1++)
    							$x1[$i1] = $x[$i1];
    						$sw = 0;
    					}
    				}
    				else
    					$ind = -1;
    			}
    			else
    				$ind = -1;
    		}
    	}
    
    	return $ind;
    }
    
    /************************/
    /* 関数値(f(x))の計算 */
    /************************/
    function snx($x, &$y)
    {
    	$y[0] = (0.5 - $x[0]) * (0.5 - $x[0]) + (1.0 - $x[1]) * (1.0 - $x[1]) - $x[2] * $x[2];
    	$y[1] = (0.0 - $x[0]) * (0.0 - $x[0]) + (1.5 - $x[1]) * (1.5 - $x[1]) - $x[2] * $x[2];
    	$y[2] = (0.5 - $x[0]) * (0.5 - $x[0]) + (2.0 - $x[1]) * (2.0 - $x[1]) - $x[2] * $x[2];
    }
    
    /*****************************************/
    /* 関数の微分の計算                      */
    /*      x : 変数値                       */
    /*      w : 微分の逆行列(後半)         */
    /*      eps : 逆行列の存在を判定する規準 */
    /*      return : =0 : 正常               */
    /*               =1 : 逆行列が存在しない */
    /*****************************************/
    function dsnx($x, &$w, $eps)
    {
    	for ($i1 = 0; $i1 < 3; $i1++) {
    		for ($i2 = 0; $i2 < 3; $i2++)
    			$w[$i1][$i2+3] = 0.0;
    		$w[$i1][$i1+3] = 1.0;
    	}
    
    	$w[0][0] = -2.0 * (0.5 - $x[0]);
    	$w[0][1] = -2.0 * (1.0 - $x[1]);
    	$w[0][2] = -2.0 * $x[2];
    	$w[1][0] = -2.0 * (0.0 - $x[0]);
    	$w[1][1] = -2.0 * (1.5 - $x[1]);
    	$w[1][2] = -2.0 * $x[2];
    	$w[2][0] = -2.0 * (0.5 - $x[0]);
    	$w[2][1] = -2.0 * (2.0 - $x[1]);
    	$w[2][2] = -2.0 * $x[2];
    
    	$sw = gauss($w, 3, 3, $eps);
    
    	return $sw;
    }
    
    /******************************************/
    /* 線形連立方程式を解く(逆行列を求める) */
    /*      $w : 方程式の左辺及び右辺         */
    /*      $n : 方程式の数                   */
    /*      $m : 方程式の右辺の列の数         */
    /*      $eps : 逆行列の存在を判定する規準 */
    /*      return : =0 : 正常                */
    /*               =1 : 逆行列が存在しない  */
    /*******************************************/
    function gauss(&$w, $n, $m, $eps)
    {
    	$ind = 0;
    	$nm  = $n + $m;
    
    	for ($i1 = 0; $i1 < $n && $ind == 0; $i1++) {
    
    		$y1 = 0.0;
    		$m1 = $i1 + 1;
    		$m2 = 0;
    				// ピボット要素の選択
    		for ($i2 = $i1; $i2 < $n; $i2++) {
    			$y2 = abs($w[$i2][$i1]);
    			if ($y1 < $y2) {
    				$y1 = $y2;
    				$m2 = $i2;
    			}
    		}
    				// 逆行列が存在しない
    		if ($y1 < $eps)
    			$ind = 1;
    				// 逆行列が存在する
    		else {
    					// 行の入れ替え
    			for ($i2 = $i1; $i2 < $nm; $i2++) {
    				$y1          = $w[$i1][$i2];
    				$w[$i1][$i2] = $w[$m2][$i2];
    				$w[$m2][$i2] = $y1;
    			}
    					// 掃き出し操作
    			$y1 = 1.0 / $w[$i1][$i1];
    
    			for ($i2 = $m1; $i2 < $nm; $i2++)
    				$w[$i1][$i2] *= $y1;
    
    			for ($i2 = 0; $i2 < $n; $i2++) {
    				if ($i2 != $i1) {
    					for ($i3 = $m1; $i3 < $nm; $i3++)
    						$w[$i2][$i3] -= $w[$i2][$i1] * $w[$i1][$i3];
    				}
    			}
    		}
    	}
    
    	return $ind;
    }
    
    ?>
    			

  5. Ruby

    ############################################
    # ニュートン法
    # (exp(x)-3.0*x=0の根)
    #      coded by Y.Suganuma
    ############################################
    
    ############################################
    # Newton法による非線形方程式(f(x)=0)の解
    #      x0 : 初期値
    #      eps1 : 終了条件1(|x(k+1)-x(k)|<eps1)
    #      eps2 : 終了条件2(|f(x(k))|<eps2)
    #      max : 最大試行回数
    #      ind : 実際の試行回数
    #            (負の時は解を得ることができなかった)
    #      fn : f(x)とその微分を計算する関数名
    #      return : 解
    #      coded by Y.Suganuma
    ############################################
    def newton(x0, eps1, eps2, max, ind, &fn) 
    
    	x1     = x0
    	x      = x1
    	ind[0] = 0
    	sw     = 0
    
    	while sw == 0 and ind[0] >= 0 
    
    		sw      = 1
    		ind[0] += 1
    		g       = fn.call(0, x1)
    
    		if g.abs() > eps2 
    			if ind[0] <= max 
    				dg = fn.call(1, x1)
    				if dg.abs() > eps2 
    					x = x1 - g / dg
    					if (x-x1).abs() > eps1 && (x-x1).abs() > eps1*x.abs() 
    						x1 = x
    						sw = 0
    					end
    				else 
    					ind[0] = -1
    				end
    			else 
    				ind[0] = -1
    			end
    		end
    	end
    
    	return x
    end
    
    #################################
    # 関数値(f(x))とその微分の計算
    #################################
    snx = Proc.new { |sw, x|
    	if sw == 0
    		Math.exp(x) - 3.0 * x
    	else
    		Math.exp(x) - 3.0
    	end
    }
    
    			# データの設定
    ind  = [0]
    eps1 = 1.0e-7
    eps2 = 1.0e-10
    max  = 20
    x0   = 0.0
    			# 実行と結果
    x = newton(x0, eps1, eps2, max, ind, &snx)
    
    print("ind=", ind[0], "  x=", x, "  f=", snx.call(0,x), "  df=", snx.call(1,x), "\n")
    			
    多次元の場合
    ############################################
    # 多次元ニュートン法
    # (3点(0.5,1.0),(0.0,1.5),(0.5,2.0)
    #   を通る円の中心と半径)
    #      coded by Y.Suganuma
    ############################################
    
    ############################################
    # Newton法による非線形方程式(f(x)=0)の解
    #      x0 : 初期値
    #      eps1 : 終了条件1(|x(k+1)-x(k)|<eps1)
    #      eps2 : 終了条件2(|f(x(k))|<eps2)
    #      max : 最大試行回数
    #      w : f(x)の微分の逆行列を入れる (n x 2n) 
    #      n : xの次数
    #      x : 解
    #      fn : f(x)とその微分を計算する関数名
    #      return : 実際の試行回数
    #            (負の時は解を得ることができなかった)
    #      coded by Y.Suganuma
    ############################################
    
    def newton(x0, eps1, eps2, max, w, n, x, &fn) 
    
    	ind = 0
    	sw  = 0
    	g   = Array.new()
    	x1  = Array.new()
    	for i1 in 0 ... n 
    		x1[i1] = x0[i1]
    		x[i1]  = x1[i1]
    	end
    
    	while sw == 0 && ind >= 0 
    
    		sw   = 1
    		ind += 1
    		fn.call(0, x1, g, w, eps2)
    
    		sw1 = 0
    		for i1 in 0 ... n 
    			if g[i1].abs() > eps2 
    				sw1 = 1
    				break
    			end
    		end
    		if sw1 > 0 
    			if ind <= max 
    				sw1 = fn.call(1, x1, g, w, eps2)
    				if sw1 == 0 
    					for i1 in 0 ... n 
    						x[i1] = x1[i1]
    						for i2 in 0 ... n 
    							x[i1] -= w[i1][i2+n] * g[i2]
    						end
    					end
    					sw1 = 0
    					for i1 in 0 ... n 
    						if (x[i1]-x1[i1]).abs() > eps1 && (x[i1]-x1[i1]).abs() > eps1*x[i1].abs() 
    							sw1 = 1
    							break
    						end
    					end
    					if sw1 > 0 
    						for i1 in 0 ... n 
    							x1[i1] = x[i1]
    						end
    						sw = 0
    					end
    				else 
    					ind = -1
    				end
    			else 
    				ind = -1
    			end
    		end
    	end
    
    	return ind
    end
    
    ############################################
    # 線形連立方程式を解く(逆行列を求める)
    #      w : 方程式の左辺及び右辺
    #      n : 方程式の数
    #      m : 方程式の右辺の列の数
    #      eps : 逆行列の存在を判定する規準
    #      return : =0 : 正常
    #               =1 : 逆行列が存在しない
    #      coded by Y.Suganuma
    ############################################
    
    def gauss(w, n, m, eps) 
    
    	nm  = n + m
    	ind = 0
    
    	for i1 in 0 ... n 
    
    		y1 = 0.0
    		m1 = i1 + 1
    		m2 = 0
    					# ピボット要素の選択
    		for i2 in i1 ... n 
    			y2 = w[i2][i1].abs()
    			if y1 < y2 
    				y1 = y2
    				m2 = i2
    			end
    		end
    					# 逆行列が存在しない
    		if y1 < eps 
    			ind = 1
    			break
    					# 逆行列が存在する
    		else    # 行の入れ替え
    			for i2 in i1 ... nm 
    				y1        = w[i1][i2]
    				w[i1][i2] = w[m2][i2]
    				w[m2][i2] = y1
    			end
    						# 掃き出し操作
    			y1 = 1.0 / w[i1][i1]
    
    			for i2 in m1 ... nm 
    				w[i1][i2] *= y1
    			end
    
    			for i2 in 0 ... n 
    				if i2 != i1 
    					for i3 in m1 ... nm 
    						w[i2][i3] -= (w[i2][i1] * w[i1][i3])
    					end
    				end
    			end
    		end
    	end
    
    	return ind
    end
    
    ########################################
    # 関数値(f(x))とその微分の計算
    #      x : 変数値
    #      y : 関数値
    #      w : 微分の逆行列(後半)
    #      eps : 逆行列の存在を判定する規準
    #      return : =0 : 正常
    #               =1 : 逆行列が存在しない
    ########################################
    snx = Proc.new { |sw, x, y, w, eps| 
    	if sw == 0   # 関数値
    		y[0] = (0.5 - x[0]) * (0.5 - x[0]) + (1.0 - x[1]) * (1.0 - x[1]) - x[2] * x[2];
    		y[1] = (0.0 - x[0]) * (0.0 - x[0]) + (1.5 - x[1]) * (1.5 - x[1]) - x[2] * x[2];
    		y[2] = (0.5 - x[0]) * (0.5 - x[0]) + (2.0 - x[1]) * (2.0 - x[1]) - x[2] * x[2];
    	else   # 微分
    		for i1 in 0 ... 3 
    			for i2 in 0 ... 3 
    				w[i1][i2+3] = 0.0
    			end
    			w[i1][i1+3] = 1.0
    		end
    
    		w[0][0] = -2.0 * (0.5 - x[0])
    		w[0][1] = -2.0 * (1.0 - x[1])
    		w[0][2] = -2.0 * x[2]
    		w[1][0] = -2.0 * (0.0 - x[0])
    		w[1][1] = -2.0 * (1.5 - x[1])
    		w[1][2] = -2.0 * x[2]
    		w[2][0] = -2.0 * (0.5 - x[0])
    		w[2][1] = -2.0 * (2.0 - x[1])
    		w[2][2] = -2.0 * x[2]
    		sw = gauss(w, 3, 3, eps)
    	end
    }
    
    			# データの設定
    eps1 = 1.0e-7
    eps2 = 1.0e-10
    max  = 20
    n    = 3
    x    = [0, 0, 0]
    x0   = [0, 0, 2]
    w    = [[0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0]]
    			# 実行と結果
    ind  = newton(x0, eps1, eps2, max, w, n, x, &snx)
    
    print("ind = ", ind, "\n")
    print("x ", x[0], " ", x[1], " ", x[2], "\n")
    snx.call(0, x, x0, w, eps1)
    print("f ", x0[0], " ", x0[1], " ", x0[2], "\n")
    			

  6. Python

    # -*- coding: UTF-8 -*-
    import numpy as np
    from math import *
    
    ############################################
    # Newton法による非線形方程式(f(x)=0)の解
    #      fn : f(x)を計算する関数名
    #      dfn : f(x)の微分を計算する関数名
    #      x0 : 初期値
    #      eps1 : 終了条件1(|x(k+1)-x(k)|<eps1)
    #      eps2 : 終了条件2(|f(x(k))|<eps2)
    #      max : 最大試行回数
    #      ind : 実際の試行回数
    #            (負の時は解を得ることができなかった)
    #      return : 解
    #      coded by Y.Suganuma
    ############################################
    
    def newton(fn, dfn, x0, eps1, eps2, max, ind) :
    
    	x1     = x0
    	x      = x1
    	ind[0] = 0
    	sw     = 0
    
    	while sw == 0 and ind[0] >= 0 :
    
    		sw      = 1
    		ind[0] += 1
    		g       = fn(x1)
    
    		if abs(g) > eps2 :
    			if ind[0] <= max :
    				dg = dfn(x1)
    				if abs(dg) > eps2 :
    					x = x1 - g / dg
    					if abs(x-x1) > eps1 and abs(x-x1) > eps1*abs(x) :
    						x1 = x
    						sw = 0
    				else :
    					ind[0] = -1
    			else :
    				ind[0] = -1
    
    	return x
    
    #######################
    # 関数値(f(x))の計算
    #######################
    def snx(x) :
    	return exp(x) - 3.0 * x
    
    ###################
    # 関数の微分の計算
    ###################
    def dsnx(x) :
    	return exp(x) - 3.0
    
    ############################################
    # ニュートン法
    # (exp(x)-3.0*x=0の根)
    #      coded by Y.Suganuma
    ############################################
    
    			# データの設定
    ind  = [0]
    eps1 = 1.0e-7
    eps2 = 1.0e-10
    max  = 20
    x0   = 0.0
    			# 実行と結果
    x = newton(snx, dsnx, x0, eps1, eps2, max, ind)
    
    print("ind=" + str(ind[0]) + "  x=" + str(x) + "  f=" + str(snx(x)) + "  df=" + str(dsnx(x)))
    			
    多次元の場合
    # -*- coding: UTF-8 -*-
    import numpy as np
    from math import *
    
    ############################################
    # Newton法による非線形方程式(f(x)=0)の解
    #      fn : f(x)を計算する関数名
    #      dfn : f(x)の微分を計算する関数名
    #      x0 : 初期値
    #      eps1 : 終了条件1(|x(k+1)-x(k)|<eps1)
    #      eps2 : 終了条件2(|f(x(k))|<eps2)
    #      max : 最大試行回数
    #      w : f(x)の微分の逆行列を入れる (n x 2n) 
    #      n : xの次数
    #      x : 解
    #      return : 実際の試行回数
    #            (負の時は解を得ることができなかった)
    #      coded by Y.Suganuma
    ############################################
    
    def newton(fn, dfn, x0, eps1, eps2, max, w, n, x) :
    
    	ind = 0
    	sw  = 0
    	g   = np.empty(n, np.float)
    	x1  = np.empty(n, np.float)
    	for i1 in range(0, n) :
    		x1[i1] = x0[i1]
    		x[i1]  = x1[i1]
    
    	while sw == 0 and ind >= 0 :
    
    		sw   = 1
    		ind += 1
    		fn(x1, g)
    
    		sw1 = 0
    		for i1 in range(0, n) :
    			if abs(g[i1]) > eps2 :
    				sw1 = 1
    				break
    		if sw1 > 0 :
    			if ind <= max :
    				sw1 = dfn(x1, w, eps2)
    				if sw1 == 0 :
    					for i1 in range(0, n) :
    						x[i1] = x1[i1]
    						for i2 in range(0, n) :
    							x[i1] -= w[i1][i2+n] * g[i2]
    					sw1 = 0
    					for i1 in range(0, n) :
    						if abs(x[i1]-x1[i1]) > eps1 and abs(x[i1]-x1[i1]) > eps1*abs(x[i1]) :
    							sw1 = 1
    							break
    					if sw1 > 0 :
    						for i1 in range(0, n) :
    							x1[i1] = x[i1]
    						sw = 0
    				else :
    					ind = -1
    			else :
    				ind = -1
    
    	return ind
    
    ############################################
    # 線形連立方程式を解く(逆行列を求める)
    #      w : 方程式の左辺及び右辺
    #      n : 方程式の数
    #      m : 方程式の右辺の列の数
    #      eps : 逆行列の存在を判定する規準
    #      return : =0 : 正常
    #               =1 : 逆行列が存在しない
    #      coded by Y.Suganuma
    ############################################
    
    def gauss(w, n, m, eps) :
    
    	nm  = n + m
    	ind = 0
    
    	for i1 in range(0, n) :
    
    		y1 = 0.0
    		m1 = i1 + 1
    		m2 = 0
    					# ピボット要素の選択
    		for i2 in range(i1, n) :
    			y2 = abs(w[i2][i1])
    			if y1 < y2 :
    				y1 = y2
    				m2 = i2
    					# 逆行列が存在しない
    		if y1 < eps :
    			ind = 1
    			break
    					# 逆行列が存在する
    		else :   # 行の入れ替え
    			for i2 in range(i1, nm) :
    				y1        = w[i1][i2]
    				w[i1][i2] = w[m2][i2]
    				w[m2][i2] = y1
    						# 掃き出し操作
    			y1 = 1.0 / w[i1][i1]
    
    			for i2 in range(m1, nm) :
    				w[i1][i2] *= y1
    
    			for i2 in range(0, n) :
    				if i2 != i1 :
    					for i3 in range(m1, nm) :
    						w[i2][i3] -= (w[i2][i1] * w[i1][i3])
    
    	return ind
    
    #######################
    # 関数値(f(x))の計算
    #######################
    def snx(x, y) :
    	y[0] = (0.5 - x[0]) * (0.5 - x[0]) + (1.0 - x[1]) * (1.0 - x[1]) - x[2] * x[2];
    	y[1] = (0.0 - x[0]) * (0.0 - x[0]) + (1.5 - x[1]) * (1.5 - x[1]) - x[2] * x[2];
    	y[2] = (0.5 - x[0]) * (0.5 - x[0]) + (2.0 - x[1]) * (2.0 - x[1]) - x[2] * x[2];
    
    ########################################
    # 関数の微分の計算
    #      x : 変数値
    #      w : 微分の逆行列(後半)
    #      eps : 逆行列の存在を判定する規準
    #      return : =0 : 正常
    #               =1 : 逆行列が存在しない
    ########################################
    def dsnx(x, w, eps) :
    	for i1 in range(0, 3) :
    		for i2 in range(0, 3) :
    			w[i1][i2+3] = 0.0
    		w[i1][i1+3] = 1.0
    
    	w[0][0] = -2.0 * (0.5 - x[0])
    	w[0][1] = -2.0 * (1.0 - x[1])
    	w[0][2] = -2.0 * x[2]
    	w[1][0] = -2.0 * (0.0 - x[0])
    	w[1][1] = -2.0 * (1.5 - x[1])
    	w[1][2] = -2.0 * x[2]
    	w[2][0] = -2.0 * (0.5 - x[0])
    	w[2][1] = -2.0 * (2.0 - x[1])
    	w[2][2] = -2.0 * x[2]
    
    	sw = gauss(w, 3, 3, eps)
    
    	return sw
    
    ############################################
    # 多次元ニュートン法
    # (3点(0.5,1.0),(0.0,1.5),(0.5,2.0)
    #   を通る円の中心と半径)
    #      coded by Y.Suganuma
    ############################################
    
    			# データの設定
    eps1 = 1.0e-7
    eps2 = 1.0e-10
    max  = 20
    n    = 3
    x    = np.array([0, 0, 0], np.float)
    x0   = np.array([0, 0, 2], np.float)
    w    = np.empty((n, 2*n), np.float)
    			# 実行と結果
    ind = newton(snx, dsnx, x0, eps1, eps2, max, w, n, x)
    
    print("ind = " + str(ind))
    print("x " + str(x[0]) + " " + str(x[1]) + " " + str(x[2]))
    snx(x, x0)
    print("f " + str(x0[0]) + " " + str(x0[1]) + " " + str(x0[2]))
    			

  7. C#

    /****************************/
    /* ニュートン法             */
    /* (exp(x)-3.0*x=0の根)   */
    /*      coded by Y.Suganuma */
    /****************************/
    using System;
    
    class Program
    {
    	static void Main()
    	{
    		Test1 ts = new Test1();
    	}
    }
    
    class Test1
    {
    	public Test1()
    	{
    		double eps1 = 1.0e-7;
    		double eps2 = 1.0e-10;
    		double x0   = 0.0;
    		int max     = 30;
    		int ind     = 0;
    					// ニュートン法の実行
    		double x = newton(x0, eps1, eps2, max, ref ind, snx, dsnx);
    					// 出力
    		Console.WriteLine("ind=" + ind + "  x=" + x + "  f=" + snx(x) + "  df=" + dsnx(x));
    	}
    					// 関数値(f(x))の計算
    	double snx(double x)
    	{
    		return Math.Exp(x) - 3.0 * x;
    	}
    					// 関数の微分の計算
    	double dsnx(double x)
    	{
    		return Math.Exp(x) - 3.0;
    	}
    
    	/*****************************************************/
    	/* Newton法による非線形方程式(f(x)=0)の解            */
    	/*      x1 : 初期値                                  */
    	/*      eps1 : 終了条件1(|x(k+1)-x(k)|<eps1)   */
    	/*      eps2 : 終了条件2(|f(x(k))|<eps2)       */
    	/*      max : 最大試行回数                           */
    	/*      ind : 実際の試行回数                         */
    	/*            (負の時は解を得ることができなかった) */
    	/*      fn : 関数値を計算する関数                    */
    	/*      dfn : 関数の微分値を計算する関数             */
    	/*      return : 解                                  */
    	/*****************************************************/
    	double newton(double x1, double eps1, double eps2, int max, ref int ind,
    	              Func<double, double> fn, Func<double, double> dfn)
    	{
    		double x = x1;
    		int sw   = 0;
    		ind      = 0;
    
    		while (sw == 0 && ind >= 0) {
    
    			ind++;
    			sw        = 1;
    			double g  = fn(x1);
    
    			if (Math.Abs(g) > eps2) {
    				if (ind <= max) {
    					double dg = dfn(x1);
    					if (Math.Abs(dg) > eps2) {
    						x = x1 - g / dg;
    						if (Math.Abs(x-x1) > eps1 && Math.Abs(x-x1) > eps1*Math.Abs(x)) {
    							x1 = x;
    							sw = 0;
    						}
    					}
    					else
    						ind = -1;
    				}
    				else
    					ind = -1;
    			}
    		}
    
    		return x;
    	}
    }
    			
    多次元の場合
    /***************************************/
    /* 多次元ニュートン法                  */
    /* (3点(0.5,1.0),(0.0,1.5),(0.5,2.0) */
    /*   を通る円の中心と半径)            */
    /*      coded by Y.Suganuma            */
    /***************************************/
    using System;
    
    class Program
    {
    	static void Main()
    	{
    		Test1 ts = new Test1();
    	}
    }
    
    class Test1
    {
    	public Test1()
    	{
    		double[] x   = new double [3];
    		double[] x0  = {0.0, 0.0, 2.0};
    		double[][] w = new double [3][];
    		w[0] = new double [6];
    		w[1] = new double [6];
    		w[2] = new double [6];
    		double eps1 = 1.0e-7;
    		double eps2 = 1.0e-10;
    		int max     = 20;
    
    		int ind = newton_m(x0, eps1, eps2, max, w, 3, x, snx, dsnx);
    
    		Console.WriteLine("ind = " + ind);
    		Console.Write("x");
    		for (int i1 = 0; i1 < 3; i1++)
    			Console.Write(" " + x[i1]);
    		Console.WriteLine();
    		Console.Write("f(x)");
    		snx(x, x0);
    		for (int i1 = 0; i1 < 3; i1++)
    			Console.Write(" " + x0[i1]);
    		Console.WriteLine();
    	}
    					// 関数値(f(x))の計算
    	int snx(double[] x, double[] y)
    	{
    		y[0] = (0.5 - x[0]) * (0.5 - x[0]) + (1.0 - x[1]) * (1.0 - x[1]) - x[2] * x[2];
    		y[1] = (0.0 - x[0]) * (0.0 - x[0]) + (1.5 - x[1]) * (1.5 - x[1]) - x[2] * x[2];
    		y[2] = (0.5 - x[0]) * (0.5 - x[0]) + (2.0 - x[1]) * (2.0 - x[1]) - x[2] * x[2];
    		return 0;
    	}
    					// 関数の微分の計算
    	int dsnx(double[] x, double[][] w, double eps)
    	{
    		for (int i1 = 0; i1 < 3; i1++) {
    			for (int i2 = 0; i2 < 3; i2++)
    				w[i1][i2+3] = 0.0;
    			w[i1][i1+3] = 1.0;
    		}
    
    		w[0][0] = -2.0 * (0.5 - x[0]);
    		w[0][1] = -2.0 * (1.0 - x[1]);
    		w[0][2] = -2.0 * x[2];
    		w[1][0] = -2.0 * (0.0 - x[0]);
    		w[1][1] = -2.0 * (1.5 - x[1]);
    		w[1][2] = -2.0 * x[2];
    		w[2][0] = -2.0 * (0.5 - x[0]);
    		w[2][1] = -2.0 * (2.0 - x[1]);
    		w[2][2] = -2.0 * x[2];
    
    		int ind = gauss(w, 3, 3, eps);
    
    		return ind;
    	}
    
    	/*****************************************************/
    	/* Newton法による多次元非線形方程式(f(x)=0)の解      */
    	/*      x1 : 初期値                                  */
    	/*      eps1 : 終了条件1(|x(k+1)-x(k)|<eps1)   */
    	/*      eps2 : 終了条件2(|f(x(k))|<eps2)       */
    	/*      max : 最大試行回数                           */
    	/*      kn1 : 関数を計算するクラスオブジェクト       */
    	/*      kn2 : 関数の微分を計算するクラスオブジェクト */
    	/*      w : f(x)の微分の逆行列を入れる (n x 2n)      */
    	/*      n : xの次数                                  */
    	/*      x : 解                                       */
    	/*      return : 実際の試行回数                      */
    	/*            (負の時は解を得ることができなかった) */
    	/*****************************************************/
    	int newton_m(double[] x0, double eps1, double eps2, int max, double[][] w,
    	             int n, double[] x,
    	             Func<double[], double[], int> fn,
    	             Func<double[], double[][], double, int> dfn)
    	{
    		double[] g  = new double [n];
    		double[] x1 = new double [n];
    		for (int i1 = 0; i1 < n; i1++) {
    			x1[i1] = x0[i1];
    			x[i1]  = x1[i1];
    		}
    		int ind = 0;
    		int sw  = 0;
    
    		while (sw == 0 && ind >= 0) {
    
    			sw = 1;
    			ind++;
    			fn(x1, g);
    
    			int sw1 = 0;
    			for (int i1 = 0; i1 < n && sw1 == 0; i1++) {
    				if (Math.Abs(g[i1]) > eps2)
    					sw1 = 1;
    			}
    			if (sw1 > 0) {
    				if (ind <= max) {
    					sw1 = dfn(x1, w, eps2);
    					if (sw1 == 0) {
    						for (int i1 = 0; i1 < n; i1++) {
    							x[i1] = x1[i1];
    							for (int i2 = 0; i2 < n; i2++)
    								x[i1] -= w[i1][i2+n] * g[i2];
    						}
    						sw1 = 0;
    						for (int i1 = 0; i1 < n && sw1 == 0; i1++) {
    							if (Math.Abs(x[i1]-x1[i1]) > eps1 &&
                                    Math.Abs(x[i1]-x1[i1]) > eps1*Math.Abs(x[i1]))
    								sw1 = 1;
    						}
    						if (sw1 > 0) {
    							for (int i1 = 0; i1 < n; i1++)
    								x1[i1] = x[i1];
    							sw = 0;
    						}
    					}
    					else
    						ind = -1;
    				}
    				else
    					ind = -1;
    			}
    		}
    
    		return ind;
    	}
    
    	/*******************************************************/
    	/* 線形連立方程式を解く(逆行列を求める)              */
    	/*      w : 方程式の左辺及び右辺                       */
    	/*      n : 方程式の数                                 */
    	/*      m : 方程式の右辺の列の数                       */
    	/*      eps : 逆行列の存在を判定する規準               */
    	/*      return : =0 : 正常                             */
    	/*               =1 : 逆行列が存在しない               */
    	/*******************************************************/
    	int gauss(double[][] w, int n, int m, double eps)
    	{
    		int ind = 0;
    		int nm  = n + m;
    
    		for (int i1 = 0; i1 < n && ind == 0; i1++) {
    
    			double y1 = .0;
    			int m1    = i1 + 1;
    			int m2    = 0;
    						// ピボット要素の選択
    			for (int i2 = i1; i2 < n; i2++) {
    				double y2 = Math.Abs(w[i2][i1]);
    				if (y1 < y2) {
    					y1 = y2;
    					m2 = i2;
    				}
    			}
    						// 逆行列が存在しない
    			if (y1 < eps)
    				ind = 1;
    						// 逆行列が存在する
    			else {
    							// 行の入れ替え
    				for (int i2 = i1; i2 < nm; i2++) {
    					y1        = w[i1][i2];
    					w[i1][i2] = w[m2][i2];
    					w[m2][i2] = y1;
    				}
    							// 掃き出し操作
    				y1 = 1.0 / w[i1][i1];
    
    				for (int i2 = m1; i2 < nm; i2++)
    					w[i1][i2] *= y1;
    
    				for (int i2 = 0; i2 < n; i2++) {
    					if (i2 != i1) {
    						for (int i3 = m1; i3 < nm; i3++)
    							w[i2][i3] -= w[i2][i1] * w[i1][i3];
    					}
    				}
    			}
    		}
    
    		return ind;
    	}
    }
    			

  8. VB

    ''''''''''''''''''''''''''''
    ' ニュートン法             '
    ' (exp(x)-3.0*x=0の根)   '
    '      coded by Y.Suganuma '
    ''''''''''''''''''''''''''''
    Module Test
    
    	Sub Main()
    
    		Dim eps1 As Double = 1.0e-7
    		Dim eps2 As Double = 1.0e-10
    		Dim x0 As Double   = 0.0
    		Dim max As Integer = 30
    		Dim ind As Integer = 0
    					' 関数値(f(x))の計算(ラムダ式)
    		Dim snx = Function(v) As Double
    		              Return Math.Exp(v) - 3.0 * v
    		          End Function
    					' 関数の微分の計算(ラムダ式)
    		Dim dsnx = Function(v) As Double
    		               Return Math.Exp(v) - 3.0
    		           End Function
    					' ニュートン法の実行
    		Dim x As Double = newton(x0, eps1, eps2, max, ind, snx, dsnx)
    					' 出力
    		Console.WriteLine("ind=" & ind & "  x=" & x & "  f=" & snx(x) & "  df=" & dsnx(x))
    	End Sub
    
    	'''''''''''''''''''''''''''''''''''''''''''''''''''''
    	' Newton法による非線形方程式(f(x)=0)の解            '
    	'      x1 : 初期値                                  '
    	'      eps1 : 終了条件1(|x(k+1)-x(k)|<eps1)   '
    	'      eps2 : 終了条件2(|f(x(k))|<eps2)       '
    	'      max : 最大試行回数                           '
    	'      ind : 実際の試行回数                         '
    	'            (負の時は解を得ることができなかった) '
    	'      fn : 関数値を計算する関数                    '
    	'      dfn : 関数の微分値を計算する関数             '
    	'      return : 解                                  '
    	'''''''''''''''''''''''''''''''''''''''''''''''''''''
    	Function newton(x1 As Double, eps1 As Double, eps2 As Double, max As Integer,
    	              ByRef ind As Integer,
    	              fn As Func(Of Double, Double), dfn As Func(Of Double, Double))
    
    		Dim x As Double   = x1
    		Dim sw As Integer = 0
    		ind = 0
    
    		Do While sw = 0 and ind >= 0
    
    			ind += 1
    			sw   = 1
    			Dim g As Double = fn(x1)
    
    			If Math.Abs(g) > eps2
    				If ind <= max
    					Dim dg As Double = dfn(x1)
    					if Math.Abs(dg) > eps2
    						x = x1 - g / dg
    						If Math.Abs(x-x1) > eps1 and Math.Abs(x-x1) > eps1*Math.Abs(x)
    							x1 = x
    							sw = 0
    						End If
    					Else
    						ind = -1
    					End If
    				Else
    					ind = -1
    				End If
    			End If
    		Loop
    
    		Return x
    
    	End Function
    
    End Module
    			
    多次元の場合
    '''''''''''''''''''''''''''''''''''''''
    ' 多次元ニュートン法                  '
    ' (3点(0.5,1.0),(0.0,1.5),(0.5,2.0) '
    '   を通る円の中心と半径)            '
    '      coded by Y.Suganuma            '
    '''''''''''''''''''''''''''''''''''''''
    Module Test
    
    	Sub Main()
    
    		Dim x(3) As Double
    		Dim x0() As Double = {0.0, 0.0, 2.0}
    		Dim w(3,6) As Double
    		Dim eps1 As Double = 1.0e-7
    		Dim eps2 As Double = 1.0e-10
    		Dim max As Integer = 20
    
    					' 関数値(f(x))の計算
    		Dim snx = Function (v1, v2)   ' 関数値の計算(ラムダ式)
    			v2(0) = (0.5 - v1(0)) * (0.5 - v1(0)) + (1.0 - v1(1)) * (1.0 - v1(1)) - v1(2) * v1(2)
    			v2(1) = (0.0 - v1(0)) * (0.0 - v1(0)) + (1.5 - v1(1)) * (1.5 - v1(1)) - v1(2) * v1(2)
    			v2(2) = (0.5 - v1(0)) * (0.5 - v1(0)) + (2.0 - v1(1)) * (2.0 - v1(1)) - v1(2) * v1(2)
    			Return 0
    		End Function
    					' 関数の微分の計算
    		Dim dsnx = Function(v1, v2, eps)   ' 関数の微分値の計算(ラムダ式)
    			For i1 As integer = 0 To 2
    				For i2 As integer = 0 To 2
    					v2(i1,i2+3) = 0.0
    				Next
    				v2(i1,i1+3) = 1.0
    			Next
    
    			v2(0,0) = -2.0 * (0.5 - v1(0))
    			v2(0,1) = -2.0 * (1.0 - v1(1))
    			v2(0,2) = -2.0 * v1(2)
    			v2(1,0) = -2.0 * (0.0 - v1(0))
    			v2(1,1) = -2.0 * (1.5 - v1(1))
    			v2(1,2) = -2.0 * v1(2)
    			v2(2,0) = -2.0 * (0.5 - v1(0))
    			v2(2,1) = -2.0 * (2.0 - v1(1))
    			v2(2,2) = -2.0 * v1(2)
    
    			Dim ind1 As Integer = gauss(v2, 3, 3, eps)
    
    			Return ind1
    		End Function
    
    		Dim ind As integer = newton_m(x0, eps1, eps2, max, w, 3, x, snx, dsnx)
    
    		Console.WriteLine("ind = " & ind)
    		Console.Write("x")
    		For i1 As Integer = 0 To 2
    			Console.Write(" " & x(i1))
    		Next
    		Console.WriteLine()
    		Console.Write("f(x)")
    		snx(x, x0)
    		For i1 As Integer = 0 To 2
    			Console.Write(" " & x0(i1))
    		Next
    		Console.WriteLine()
    
    	End Sub
    
    	'''''''''''''''''''''''''''''''''''''''''''''''''''''
    	' Newton法による多次元非線形方程式(f(x)=0)の解      '
    	'      x1 : 初期値                                  '
    	'      eps1 : 終了条件1(|x(k+1)-x(k)|<eps1)   '
    	'      eps2 : 終了条件2(|f(x(k))|<eps2)       '
    	'      max : 最大試行回数                           '
    	'      kn1 : 関数を計算するクラスオブジェクト       '
    	'      kn2 : 関数の微分を計算するクラスオブジェクト '
    	'      w : f(x)の微分の逆行列を入れる (n x 2n)      '
    	'      n : xの次数                                  '
    	'      x : 解                                       '
    	'      return : 実際の試行回数                      '
    	'            (負の時は解を得ることができなかった) '
    	'''''''''''''''''''''''''''''''''''''''''''''''''''''
    	Function newton_m(x0() As Double, eps1 As Double, eps2 As Double, max As Integer,
    	             w(,) As Double, n As Integer, x() As Double,
    	             fn As Func(Of Double(), Double(), Integer),
    	             dfn As Func(Of Double(), Double(,), Double, Integer))
    
    		Dim g(n) As Double
    		Dim x1(n) As Double
    		For i1 As Integer = 0 To n-1
    			x1(i1) = x0(i1)
    			x(i1)  = x1(i1)
    		Next
    		Dim ind As Integer = 0
    		Dim sw As Integer  = 0
    
    		Do While sw = 0 and ind >= 0
    
    			sw   = 1
    			ind += 1
    			fn(x1, g)
    
    			Dim sw1 As Integer = 0
    			Dim i1 As Integer  = 0
    			Do While i1 < n and sw1 = 0
    				If Math.Abs(g(i1)) > eps2
    					sw1 = 1
    				End If
    				i1 += 1
    			Loop
    			If sw1 > 0
    				If ind <= max
    					sw1 = dfn(x1, w, eps2)
    					If sw1 = 0
    						For i1 = 0 To n-1
    							x(i1) = x1(i1)
    							For i2 As Integer = 0 To n-1
    								x(i1) -= w(i1,i2+n) * g(i2)
    							Next
    						Next
    						sw1 = 0
    						i1  = 0
    						Do While i1 < n and sw1 = 0
    							If Math.Abs(x(i1)-x1(i1)) > eps1 and
                                   Math.Abs(x(i1)-x1(i1)) > eps1*Math.Abs(x(i1))
    								sw1 = 1
    							End If
    							i1 += 1
    						Loop
    						If sw1 > 0
    							For i1 = 0 To n-1
    								x1(i1) = x(i1)
    							Next
    							sw = 0
    						End If
    					Else
    						ind = -1
    					End If
    				Else
    					ind = -1
    				End If
    			End If
    		Loop
    
    		Return ind
    
    	End Function
    
    	''''''''''''''''''''''''''''''''''''''''''
    	' 線形連立方程式を解く(逆行列を求める) '
    	'      w : 方程式の左辺及び右辺          '
    	'      n : 方程式の数                    '
    	'      m : 方程式の右辺の列の数          '
    	'      eps : 逆行列の存在を判定する規準  '
    	'      return : =0 : 正常                '
    	'               =1 : 逆行列が存在しない  '
    	''''''''''''''''''''''''''''''''''''''''''
    	Function gauss(w(,) As Double, n As Integer, m As Integer, eps As Double) As Integer
    
    		Dim ind As Integer = 0
    		Dim nm As Integer  = n + m
    
    		Dim i1 As Integer = 0
    		Do While i1 < n and ind = 0
    
    			Dim y1 As Double  = 0.0
    			Dim m1 As Integer = i1 + 1
    			Dim m2 As Integer = 0
    						' ピボット要素の選択
    			For i2 As Integer = i1 To n-1
    				Dim y2 As Double = Math.Abs(w(i2,i1))
    				If y1 < y2
    					y1 = y2
    					m2 = i2
    				End If
    			Next
    						' 逆行列が存在しない
    			If y1 < eps
    				ind = 1
    						' 逆行列が存在する
    			Else
    							' 行の入れ替え
    				For i2 As Integer = i1 To nm-1
    					y1       = w(i1,i2)
    					w(i1,i2) = w(m2,i2)
    					w(m2,i2) = y1
    				Next
    							' 掃き出し操作
    				y1 = 1.0 / w(i1,i1)
    
    				For i2 As Integer = m1 To nm-1
    					w(i1,i2) *= y1
    				Next
    
    				For i2 As Integer = 0 To n-1
    					If i2 <> i1
    						For i3 As Integer = m1 To nm-1
    							w(i2,i3) -= w(i2,i1) * w(i1,i3)
    						Next
    					End If
    				Next
    			End If
    			i1 += 1
    		Loop
    
    		Return ind
    
    	End Function
    
    End Module
    			

静岡理工科大学 菅沼ホーム 目次 索引