| 情報学部 | 菅沼ホーム | 目次 | 索引 |
/****************************/
/* ニュートン法 */
/* （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 : 終了条件１（｜x(k+1)-x(k)｜＜eps1） */
/* eps2 : 終了条件２（｜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;
}
/***************************************/
/* 多次元ニュートン法 */
/* （３点(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 : 終了条件１（｜x(k+1)-x(k)｜＜eps1） */
/* eps2 : 終了条件２（｜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);
}
/****************************/
/* ニュートン法 */
/* （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 : 終了条件１（｜x(k+1)-x(k)｜＜eps1） */
/* eps2 : 終了条件２（｜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;
}
}
/***************************************/
/* 多次元ニュートン法 */
/* （３点(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 : 終了条件１（｜x(k+1)-x(k)｜＜eps1） */
/* eps2 : 終了条件２（｜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);
}
}
<!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 : 終了条件１（｜x(k+1)-x(k)｜＜eps1） */
/* eps2 : 終了条件２（｜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 : 終了条件１（｜x(k+1)-x(k)｜＜eps1） */
/* eps2 : 終了条件２（｜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>
<?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 : 終了条件１（｜x(k+1)-x(k)｜＜eps1） */
/* eps2 : 終了条件２（｜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
/***************************************/
/* 多次元ニュートン法 */
/* （３点(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 : 終了条件１（｜x(k+1)-x(k)｜＜eps1） */
/* eps2 : 終了条件２（｜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;
}
?>
############################################
# ニュートン法
# （exp(x)-3.0*x=0の根）
# coded by Y.Suganuma
############################################
############################################
# Newton法による非線形方程式(f(x)=0)の解
# x0 : 初期値
# eps1 : 終了条件１（｜x(k+1)-x(k)｜＜eps1）
# eps2 : 終了条件２（｜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")
############################################
# 多次元ニュートン法
# （３点(0.5,1.0),(0.0,1.5),(0.5,2.0)
# を通る円の中心と半径）
# coded by Y.Suganuma
############################################
############################################
# Newton法による非線形方程式(f(x)=0)の解
# x0 : 初期値
# eps1 : 終了条件１（｜x(k+1)-x(k)｜＜eps1）
# eps2 : 終了条件２（｜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")
# -*- coding: UTF-8 -*-
import numpy as np
from math import *
############################################
# Newton法による非線形方程式(f(x)=0)の解
# fn : f(x)を計算する関数名
# dfn : f(x)の微分を計算する関数名
# x0 : 初期値
# eps1 : 終了条件１（｜x(k+1)-x(k)｜＜eps1）
# eps2 : 終了条件２（｜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 : 終了条件１（｜x(k+1)-x(k)｜＜eps1）
# eps2 : 終了条件２（｜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
############################################
# 多次元ニュートン法
# （３点(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]))
/****************************/
/* ニュートン法 */
/* （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 : 終了条件１（｜x(k+1)-x(k)｜＜eps1） */
/* eps2 : 終了条件２（｜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;
}
}
/***************************************/
/* 多次元ニュートン法 */
/* （３点(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 : 終了条件１（｜x(k+1)-x(k)｜＜eps1） */
/* eps2 : 終了条件２（｜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;
}
}
''''''''''''''''''''''''''''
' ニュートン法 '
' （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 : 終了条件１（｜x(k+1)-x(k)｜＜eps1） '
' eps2 : 終了条件２（｜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
'''''''''''''''''''''''''''''''''''''''
' 多次元ニュートン法 '
' （３点(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 : 終了条件１（｜x(k+1)-x(k)｜＜eps1） '
' eps2 : 終了条件２（｜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
| 情報学部 | 菅沼ホーム | 目次 | 索引 |