要求された機能を満たすニュートン法のプログラム例を示す。
(1) ニュートン法により正の平方根を求めるプログラム newton.c の改良
(1-1) f (x ) と f '(x ) を独立させたプログラム
(newton1_1.c)。
/*
newton1_1.c: f(x), f'(x) 独立バージョン
*/
#include <stdio.h> // printf, fprintf, fgets, sscanf
#include <math.h> // fabs
double a = 2;
double f(double x)
{
return x * x - a; // f(x) = x2 - a
}
double df(double x)
{
return 2 * x; // f'(x) = 2x
}
int main(void)
{
int n = 0;
double x, x0, err, eps = 1.0e-10;
char s[128];
fprintf(stderr, " a = "); fgets(s,128,stdin); sscanf(s, "%lf", &a);
while (a <= 0.0) {
fprintf(stderr, "'a' には正の数を入れてください。\n");
fprintf(stderr, " a = "); fgets(s,128,stdin); sscanf(s, "%lf", &a);
}
x = (a + 1.0) / 2.0;
printf("# n, x, err\n");
printf("%4d, % .15e\n", n, x);
do {
n++;
x0 = x;
x = x0 - f(x0) / df(x0);
err = fabs(x - x0);
printf("%4d, % .15e, % .15e\n", n, x, err);
} while (err >= eps);
printf("\n# sqrt(%e) = % .15e\n", a, x);
return 0;
}
(1-2) アルゴリズム部分を newton() として独立させたプログラム
(newton1_2.c)。
/*
newton1_2.c: アルゴリズム部分独立バージョン
*/
#include <stdio.h> // printf, fprintf, fgets, sscanf
#include <math.h> // fabs
double a = 2;
double f(double x)
{
return x * x - a; // f(x) = x2 - a
}
double df(double x)
{
return 2 * x; // f'(x) = 2x
}
double newton(double x, double eps)
/*
Solving f(x) = 0 by Newton Method
x: initial value of x
eps: error
f(x): left term of equation f(x)=0
df(x): derivative of f(x)
*/
{
int n = 0;
double x0, err;
printf("# n, x, err\n");
printf("%4d, % .15e\n", n, x);
do {
n++;
x0 = x;
x = x0 - f(x0) / df(x0);
err = fabs(x - x0);
printf("%4d, % .15e, % .15e\n", n, x, err);
} while (err >= eps);
return x;
}
int main(void)
{
double x, eps = 1.0e-10;
char s[128];
fprintf(stderr, " a = "); fgets(s,128,stdin); sscanf(s, "%lf", &a);
while (a <= 0.0) {
fprintf(stderr, "'a' には正の数を入れてください。\n");
fprintf(stderr, " a = "); fgets(s,128,stdin); sscanf(s, "%lf", &a);
}
x = (a + 1.0) / 2.0;
printf("\n# sqrt(%e) = % .15e\n", a, newton(x, eps));
return 0;
}
(1-3) newton() の引数に関数ポインタを追加したプログラム
(newton1_3.c)。
/*
newton1_3.c: 関数ポインタ使用バージョン
*/
#include <stdio.h> // printf, fprintf, fgets, sscanf
#include <math.h> // fabs
double a = 2;
double f(double x)
{
return x * x - a; // f(x) = x^2 - a
}
double df(double x)
{
return 2 * x; // f'(x) = 2x
}
double newton(double x, double (*f)(double), double (*df)(double), double eps)
/*
Solving f(x) = 0 by Newton Method
x: initial value of x
f(x): left term of equation f(x)=0
df(x): derivative of f(x)
eps: error
*/
{
int n = 0;
double x0, err;
printf("# n, x, err\n");
printf("%4d, % .15e\n", n, x);
do {
n++;
x0 = x;
x = x0 - f(x0) / df(x0);
err = fabs(x - x0);
printf("%4d, % .15e, % .15e\n", n, x, err);
} while (err >= eps);
return x;
}
int main(void)
{
double x, eps = 1.0e-10;
char s[128];
fprintf(stderr, " a = "); fgets(s,128,stdin); sscanf(s, "%lf", &a);
while (a <= 0.0) {
fprintf(stderr, "'a' には正の数を入れてください。\n");
fprintf(stderr, " a = "); fgets(s,128,stdin); sscanf(s, "%lf", &a);
}
x = (a + 1.0) / 2.0;
printf("\n# sqrt(%e) = % .15e\n", a, newton(x, f, df, eps));
return 0;
}
(1-4) a の数値をコマンド行で入力するようにしたプログラム
(newton1_4.c)。
→[付録2]
/*
newton1_4.c: コマンド行入力バージョン
*/
#include <stdio.h> // printf, fprintf, fgets, sscanf
#include <math.h> // fabs, atof
double a = 2;
double f(double x)
{
return x * x - a; // f(x) = x^2 - a
}
double df(double x)
{
return 2 * x; // f'(x) = 2x
}
double newton(double x, double (*f)(double), double (*df)(double), double eps)
/*
Solving f(x) = 0 by Newton Method
x: initial value of x
f(x): left term of equation f(x)=0
df(x): derivative of f(x)
eps: error
*/
{
int n = 0;
double x0, err;
printf("# n, x, err\n");
printf("%4d, % .15e\n", n, x);
do {
n++;
x0 = x;
x = x0 - f(x0) / df(x0);
err = fabs(x - x0);
printf("%4d, % .15e, % .15e\n", n, x, err);
} while (err >= eps);
return x;
}
int main(int argc, char *argv[])
{
double x, eps = 1.0e-10;
char s[128];
if (argc < 2) {
fprintf(stderr, "usage: %s a\n", argv[0]);
return 1;
}
a = atof(argv[1]);
x = (a + 1.0) / 2.0;
printf("\n# sqrt(%e) = % .15e\n", a, newton(x, f, df, eps));
return 0;
}
(2) 2変数のニュートン法のプログラム newton2.c の改良
(2-1) f (x ) と J (x ) をまとめたプログラム
(newton2_1.c)。
→[付録2]
/*
newton2_1.c: f(x), J(x) をまとめたバージョン
x2 + y2 = 1, y = x3 -> (x, y) = (0.826, 0.5636), (-0.826, -0.5636)
*/
#include <stdio.h> // printf, fprintf, fgets, sscanf
#include <math.h> // hypot(hypotenuse): 直角三角形の斜辺を求める関数
#define ND 3
int func(double x, double y, double f[])
{
f[1] = x * x + y * y - 1; // f1(x,y) = x2 + y2 - 1
f[2] = x * x * x - y; // f2(x,y) = x3 - y
return 0;
}
int Jcalc(double x, double y, double J[][ND])
{
J[1][1] = 2 * x; // df1(x,y)/dx = 2x
J[1][2] = 2 * y; // df1(x,y)/dy = 2y
J[2][1] = 3 * x * x; // df2(x,y)/dx = 3x2
J[2][2] = -1; // df2(x,y)/dy = -1
return 0;
}
int ucalc(double x, double y, double *ux, double *uy)
/*
u = -J-1(x)*f(x)
*/
{
double det, f[ND], J[ND][ND];
func(x, y, f); // 配列f[]に f(x,y) の各成分を格納
Jcalc(x, y, J); // 配列J[][]に J(x,y) の各成分を格納
det = J[1][1] * J[2][2] - J[1][2] * J[2][1];
*ux = -1 / det * (J[2][2] * f[1] - J[1][2] * f[2]);
*uy = 1 / det * (J[2][1] * f[1] - J[1][1] * f[2]);
return 0;
}
int newton2(double *x, double *y, double eps)
{
int n = 0;
double ux, uy, err;
printf("# n, x, y, err\n");
printf("%4d, % .15e, % .15e\n", n, *x, *y);
do {
n++;
ucalc(*x, *y, &ux, &uy);
*x += ux; *y += uy; // xk+1 = xk + u(xk)
err = hypot(ux, uy);
printf("%4d, % .15e, % .15e, % .15e\n", n, *x, *y, err);
} while (err >= eps);
return 0;
}
// ------------ Main ------------
int main(void)
{
double x = 1, y = 0.5, eps = 1e-10;
char s[128];
fprintf(stderr,"# x0 y0 = "); fgets(s,128,stdin);
sscanf(s,"%lf%lf",&x,&y);
newton2(&x, &y, eps);
printf("\n# (x, y) = (% .15e, % .15e)\n", x, y);
return 0;
}
(2-2) 変数x, y を配列x[1], x[2] に、変数ux, uy を
u[1], u[2] に割り当てて newton2() の引数の数を減らしたプログラム
(newton2_2.c)。
/*
newton2_2.c: x, y を x[]、ux, uy を u[] にまとめたバージョン
x2 + y2 = 1, y = x3 -> (x, y) = (0.826, 0.5636), (-0.826, -0.5636)
*/
#include <stdio.h> // printf, fprintf, fgets, sscanf
#include <math.h> // hypot(hypotenuse): 直角三角形の斜辺を求める関数
#define ND 3
int func(double x[], double f[])
{
f[1] = x[1] * x[1] + x[2] * x[2] - 1; // f1(x,y) = x2 + y2 - 1
f[2] = x[1] * x[1] * x[1] - x[2]; // f2(x,y) = x3 - y
return 0;
}
int Jcalc(double x[], double J[][ND])
{
J[1][1] = 2 * x[1]; // df1(x,y)/dx = 2x
J[1][2] = 2 * x[2]; // df1(x,y)/dy = 2y
J[2][1] = 3 * x[1] * x[1]; // df2(x,y)/dx = 3x2
J[2][2] = -1; // df2(x,y)/dy = -1
return 0;
}
int ucalc(double x[], double u[])
/*
u = -J-1(x)*f(x)
*/
{
double det, f[ND], J[ND][ND];
func(x, f); // 配列f[]に f(x) の各成分を格納
Jcalc(x, J); // 配列J[][]に J(x) の各成分を格納
det = J[1][1] * J[2][2] - J[1][2] * J[2][1];
u[1] = -1 / det * (J[2][2] * f[1] - J[1][2] * f[2]);
u[2] = 1 / det * (J[2][1] * f[1] - J[1][1] * f[2]);
return 0;
}
int newton2(double x[], double eps)
{
int n = 0;
double u[ND], err;
printf("# n, x, y, err\n");
printf("%4d, % .15e, % .15e\n", n, x[1], x[2]);
do {
n++;
ucalc(x, u);
x[1] += u[1]; x[2] += u[2]; // xk+1 = xk + u(xk)
err = hypot(u[1],u[2]);
printf("%4d, % .15e, % .15e, % .15e\n", n, x[1], x[2], err);
} while (err >= eps);
return 0;
}
// ------------ Main ------------
int main(void)
{
double x[ND], eps = 1.0e-10;
char s[128];
x[1] = 1; x[2] = 0.5;
fprintf(stderr,"# x0 y0 = "); fgets(s,128,stdin);
sscanf(s,"%lf%lf",&x[1],&x[2]);
newton2(x, eps);
printf("\n# (x, y) = (% .15e, % .15e)\n", x[1], x[2]);
return 0;
}
(2-3) newton2()の引数に関数ポインタを追加したプログラム
(newton2_3.c)。
/*
newton2_3.c: 関数ポインタ追加バージョン
x2 + y2 = 1, y = x3 -> (x, y) = (0.826, 0.5636), (-0.826, -0.5636)
*/
#include <stdio.h> // printf, fprintf, fgets, sscanf
#include <math.h> // hypot(hypotenuse): 直角三角形の斜辺を求める関数
#define ND 3
int func(double x[], double f[])
{
f[1] = x[1] * x[1] + x[2] * x[2] - 1; // f1(x,y) = x2 + y2 - 1
f[2] = x[1] * x[1] * x[1] - x[2]; // f2(x,y) = x3 - y
return 0;
}
int Jcalc(double x[], double J[][ND])
{
J[1][1] = 2 * x[1]; // df1(x,y)/dx = 2x
J[1][2] = 2 * x[2]; // df1(x,y)/dy = 2y
J[2][1] = 3 * x[1] * x[1]; // df2(x,y)/dx = 3x2
J[2][2] = -1; // df2(x,y)/dy = -1
return 0;
}
int ucalc(double x[], double u[], int (*func)(double*, double*),
int (*Jcalc)(double*,double J[][ND]))
/*
u = -J-1(x)*f(x)
*/
{
double det, f[ND], J[ND][ND];
func(x, f); // 配列f[]に f(x) の各成分を格納
Jcalc(x, J); // 配列J[][]に J(x) の各成分を格納
det = J[1][1] * J[2][2] - J[1][2] * J[2][1];
u[1] = -1 / det * (J[2][2] * f[1] - J[1][2] * f[2]);
u[2] = 1 / det * (J[2][1] * f[1] - J[1][1] * f[2]);
return 0;
}
int newton2(double x[], int (*func)(double*,double*),
int (*Jcalc)(double*,double J[][ND]), double eps)
{
int n = 0;
double u[ND], err;
printf("# n, x, y, err\n");
printf("%4d, % .15e, % .15e\n", n, x[1], x[2]);
do {
n++;
ucalc(x, u, func, Jcalc);
x[1] += u[1]; x[2] += u[2]; // xk+1 = xk + u(xk)
err = hypot(u[1], u[2]);
printf("%4d, % .15e, % .15e, % .15e\n",n, x[1], x[2], err);
} while (err >= eps);
return 0;
}
// ------------ Main ------------
int main(void)
{
double x[ND], eps = 1e-10;
char s[128];
x[1] = 1; x[2] = 0.5;
fprintf(stderr,"# x0 y0 = "); fgets(s,128,stdin);
sscanf(s,"%lf%lf",&x[1],&x[2]);
newton2(x, func, Jcalc, eps);
printf("\n# (x, y) = (% .15e, % .15e)\n", x[1], x[2]);
return 0;
}
○新たな実行例のための f(x), J(x)
int func2(double x[], double f[])
{
f[1] = x[1] * x[1] + x[2] * x[2] / 4 - 1; // f1(x,y) = x2 + (y/2)2 - 1
f[2] = x[1] * x[1] - x[2] * x[2] + 1; // f2(x,y) = x2 - y2 + 1
return 0;
}
int Jcalc2(double x[], double J[][ND])
{
J[1][1] = 2 * x[1]; // df1(x,y)/dx = 2x
J[1][2] = x[2] / 2; // df1(x,y)/dy = y/2
J[2][1] = 2 * x[1]; // df2(x,y)/dx = 2x
J[2][2] = -2 * x[2]; // df2(x,y)/dy = -2y
return 0;
}
これを ucalc() より前に追加し、main() にて newton2(x, func2, Jcalc2, eps) と呼び出せばよい
(newton2_4.c)。
Z:\proen1\kadai4>newton2_4
# x0 y0 = 1 0.5
# n, x, y, err
0, 1.000000000000000e+000, 5.000000000000000e-001
1, 8.000000000000000e-001, 1.850000000000000e+000, 1.364734406395618e+000
2, 7.750000000000000e-001, 1.357432432432433e+000, 4.932015902442228e-001
3, 7.745967741935483e-001, 1.268064150511886e+000, 8.936919158457984e-002
4, 7.745966692414905e-001, 1.264914984197939e+000, 3.149166315695276e-003
5, 7.745966692414833e-001, 1.264911064073426e+000, 3.920124513106996e-006
6, 7.745966692414833e-001, 1.264911064067352e+000, 6.074514284039246e-012
# (x, y) = ( 7.745966692414833e-001, 1.264911064067352e+000)
【目次】 |
【1.】 |
【2.のつづき】 |
【3.】 |
【付録1】 |
【付録2】 |
【付録3】