数値解析
2023年度前期 第4週 [5月11日]

三浦 憲二郎講義アウトライン [5月11日]•非線形方程式•復習• ２分法•非線形方程式•ニュートン法 （1変数）復習：非線形方程式 p.70
•代数方程式の解•4次方程式までは公式があり，解ける．
•5次以上の方程式に代数的解法は存在しない．
•数値的に解くしかない．
•高次の連立方程式，三角関数など変数のn乗で表せない項を含む方程式•数値的に解くしかない．

復習：2分法 p.70
•定理4.1 （中間値の定理） 関数f(x)は閉区間[a,b]で連続で f(a)≠f(b) ならば，f(a) と f(b) の間の任意の数 k に対して f(c)=k となる c (a<c<b)が存在する．
2分法•中間値の定理より，f(a) f(b) < 0 ならば，f(c) =0, a < c < b となる cが存在する．
•(1) 何らかの方法で f(a) f(b) < 0 となる閉区間 [a, b] を求める．
•(2) c1 = (a+b)/2•(2a) f(a) f(c1) < 0 のとき，解αは [a, c1] に存在する．
•(2b) f(c1) f(b) < 0 のとき，解αは [c1, b] に存在する．
(2a) では [a, b] の代わりに[a,c1], (2b)では [c1, b]とする．
• ステップ(2)に戻り，閉区間が十分に小さくなるまで繰り返す

復習：2分法のアルゴリズム p.71

復習：2分法
区間 [a, b]に解が1つしか存在しないのならば，•n回目の区間•dn < εのときに計算を打ち切るとすると，必要な計算回数nは，を満たす最小の自然数
絶対値誤差
ステップ1) 区間 [a, b] の求め方•（1） 最初の区間 [xmin, xmax] ，および 微小区間の幅 h を与える．
•（2） n=(xmax-xmin) / h として分割数を定め，x0=xminとする．
•（3） k=1,2,…,nに対して，xk=xmin+khとし，f(xk-1) f(xk) < 0 なら
[xk-1, xk] を対象区間にする．

復習：2分法：プログラム#include <stdio.h>#include <math.h>double bisection(double a,double b,double eps); /* 2分法 */double f(double x); /* 関数の定義 */int main(void){
double a,b,x,h,y1,y2,eps=pow(2.0,-30.0);
int n;
printf("初期区間[a,b]を入力してください．---> a b¥n");scanf("%lf%lf",&a,&b);printf("区間の分割数nを入力してください．---> n¥n");scanf("%d",&n);}/* 対象区間を探索しながら2分法を適用 */h=(b-a)/n; y1=f(a);for(x=a+h;x<=b;x+=h){y2=f(x);if(y1*y2<0.0){
printf("求める答えはx=%fです.¥n",bisection(x-h,x,eps));}y1=y2;}
return 0;2分法：プログラム/* 2分法 */double bisection(double a,double b,double eps){
double c;do{}

c=0.5*(a+b);if(f(a)*f(c)<0){
b=c;}
else{
a=c;}
}while(fabs(b-a) >=eps)；/* fabs()は絶対値を返す．「C言語入門」p.264 */c=0.5*(a+b);
return c;/* 関数の定義 */double f(double x){
return x*(x*x*(x*x-5.0)+4.0); /* x*x*x*x*x-5.0*x*x*x+4.0*x */}復習：プログラム：実行結果
初期区間[a,b]を入力してください．- - ->a b-3 3

区間の分割数nを入力してください．- - -> n10
求める答えはx=-2.000000です．
求める答えはx=-1.000000です．
求める答えはx=-0.000000です．
求める答えはx=1.000000です．

25

求める答えはx=2.000000です．

2015
105-3-2-10-5-10-15-20-250123

反復法と縮小写像の原理 p.74
2分法 解の存在範囲を探索しながら狭めていく． 直接探索法
ニュートン法 初期値x0から出発して解αに収束するような列{xn}を作り，xnがαに十分近づいたときに計算を打ち切る．
反復法，あるいは逐次反復法• 反復法による列{xn}g(x) : 反復関数，式(1) x : 不動点, 式(2) : 不動点反復φ(x)≠0となるように選んで，

ニュートン法 p.77
ニュートン法の原理•

解αの近傍でC 2級とする．(2次導関数が連続)•テイラー展開•

ニュートン反復列

ニュートン法 p.77
収束しない場合がある．

ニュートン法：プログラム#include <stdio.h>#include <math.h>#define EPS pow(10.0,-8.0) /* epsilonの設定 */#define NMAX 10 /* 最大反復回数 */*/void newton(double x); /* Newton法
double f(double x);/* f(x)の計算 */double df(double x);/* f'(x)の計算 */int main(void){
double x;
printf("初期値x0を入力してください¥n");scanf("%lf",&x);
newton(x);return 0;}
void newton(double x){/* Newton法 */int n=0; double d;do{d=-f(x)/df(x);x=x+d;
n++;}while(fabs(d)>EPS&&n<NMAX);if(n==NMAX){
printf("答えは見つかりませんでした¥n");}else{
printf("答えはx=%fです．¥n",x);}}/* 関数の定義 */double f(double x){
return x-cos(x);}double df(double x){
return 1.0+sin(x);}ニュートン法：実行結果 p.79
実行結果

x-cos(x)=0(1回目)
初期値x0を入力してください．
３
答えはx=0.739085です．
（2回目)
初期値x0を入力してください．
4
答えがみつかりませんでした．

6543（3回目）
初期値x0を入力してください．
５
答えがみつかりませんでした．

210-1-20

1234567

収束の速さ p.82
(１） 線形収束
αに収束する反復列{xn}がAＡはnに依存しない定数(２） p次収束

AMはnに依存しない定数収束の速さ p.82
定理4.3 f(x)=0の解ｘ＝αが単解のとき，ニュートン法は２次収束する．
（証明）
ｆ（α）＝0， ｆ’（α）≠0 に注意すると，テイラーの公式より

なぜならば，収束の速さ p.82
ニュートン法が収束するような区間において，となるような定数A, Bを選べば，が成り立つ．これはニュートン法が（収束するならば）2次収束することを示している．Excelによるグラフの作成
2次関数のグラフ#include <stdio.h>#include <stdlib.h>
int main(void){

int i; double x,y;FILE *fout;if((fout=fopen("output.dat","w"))==NULL){

printf("ファイルが作成できません:output.dat ¥n");}exit(1);for(x=-3;x<=3;x+=0.2){
y=x*x;}}

fprintf(fout,"%lf %lf¥n",x,y);

fclose(fout);Excelによるグラフの作成
output.dat-39

グラフ：散布図-2.8 7.84-2.6 6.76-2.4 5.76

10-28

9-2.2 4.84
4-1.8 3.24

7-1.6 2.56

6-1.4 1.96

5-1.2 1.44-14

1-0.8 0.64

3-0.4 0.16

1-0.6 0.36

-0.2 0.04

0…

0

系列12-4-3-2-100

1234まとめ•非線形方程式•ニュートン法•C言語の基礎•ファイル入出力•Excelによるグラフの作成