［更新日］1998.11.06
問１３ 次のFortranプログラムの説明及びプログラムを読んで、 設問１〜３に答えよ。
〔プログラムの説明〕
n 次方程式の次数，係数及び初期値を入力して，n 次方程式（1 ≦ n ≦ 10）
an xn ＋ an−1 xn−1 ＋ … ＋ a1 x ＋ a0 ＝ 0
の解をニュートン法によって求めるプログラムである。
ここでは，
(1) f (x) ＝ an xn ＋ an−1 xn−1 ＋ … ＋ a1 x ＋ a0 とする。
(2) f ' (x) ＝ n an xn−1 ＋ (n−1) an−1 xn−2 ＋ … ＋ a1 である。
(3) 初期値を x1 として，k 回繰り返して計算した結果得られる 近似値を xk＋1 とするとき，それを次式で求める。
そして，｜xk＋1 − xk ｜＜ε（ε：収束判定定数）となったとき収束したと見なし，近似値 xk＋1 を 解として出力する。
プログラムでは，εを 10−6，反復回数の最大を 10 回とし， 収束しなかった場合は“no convergence”を出力する。
計算例として，方程式
f (x) ＝ x3 −8.335837 x2 ＋ 20.74309 x − 15.49191 ＝ 0
の解を，初期値を 1.0 として求めたときの出力結果を次に示す。
k x(k+1) x(k) f(x(k))/f'(x(k)) 1 1.2948006 1.0000000 -0.2948006 2 1.3994567 1.2948006 -0.1046561 3 1.4139398 1.3994567 -0.0144831 4 1.4142132 1.4139398 -0.0002733 5 1.4142132 1.4142132 0.0000000 solution = 1.4142132
なお，この方程式の三つの解は，1.414213 …，2.449489 …，4.472135 … である。
また，24 行目と 25 行目の間に xk1，fx，dfx の値を出力する文を挿入して実行したところ，次の出力が得られた。
xk1 fx dfx 1.2948006 -2.0846577 7.0714149 1.3994567 -0.4381027 4.1861191 1.4139398 -0.0476093 3.2872391 1.4142132 -0.0008659 3.1680222 1.4142132 0.0000000 3.1657848
〔プログラム〕
（行番号） 01 real,dimension(0:10):: a 02 real:: xk,xk1,fd 03 real:: fx,dfx 04 integer:: n,k,icnt,icnv 05 real,parameter:: eps=1.0e-06 06 integer,parameter:: icmax=10 07 read(*,*) n,(a(k),k=n,0,-1) 08 read(*,*) xk 09 icnt = 0 10 icnv = 0 11 write(*,200) 12 do while((icnt<icmax).and.(icnv==0)) 13 icnt = icnt+1 14 fx = a(n) 15 dfx = real(n)*a(n) 16 do k = n-1,0,-1 17 fx = fx*xk+a(k) 18 if(k>0) then 19 dfx = dfx*xk+real(k)*a(k) 20 end if 21 end do 22 fd = fx/dfx 23 xk1 = xk-fd 24 write(*,210) icnt,xk1,xk,fd 25 if(abs(xk1-xk)<eps) then 26 icnv = 1 27 end if 28 xk = xk1 29 end do 30 if(icnv==1) then 31 write(*,220) xk 32 else 33 write(*,230) 34 end if 35 200 format(' k x(k+1) x(k) '& 36 & 'f(x(k))/f''(x(k))') 37 210 format(i5,3f16.7) 38 220 format('solution = ',f16.7) 39 230 format('no convergence') 40 end
設問１ このプログラムの 16 行目から 21 行目までを変更すれば，18 行目から 20 行目までの IF 構文を用いなくても，同じ結果を得るプログラムに書き換えられる。 この変更として正しい答えを，解答群の中から選べ。
解答群 ア do k = n-1,1,-1 イ do k = n-1,0,-1 fx = fx*xk+a(k) fx = fx*xk+a(k) dfx = dfx*xk+real(k)*a(k) dfx = dfx*xk+real(k)*a(k) end do end do fx = fx*xk+a(1) fx = fx*xk+a(1) ウ do k = n-1,1,-1 エ do k = n-1,0,-1 fx = fx*xk+a(k) fx = fx*xk+a(k) dfx = dfx*xk+real(k)*a(k) dfx = dfx*xk+real(k)*a(k) end do end do fx = fx*xk+a(0) fx = fx*xk+a(0)
設問２ このプログラムの 14 行目から 21 行目までを取り出して，fx/dfx を 計算する関数副プログラム funcfdfx を作成した。関数副プログラム funcfdfx を 利用するために，このプログラムの 3 行目及び 14 行目から 21 行目までを 削除したとき，22 行目の変更として正しい答えを，解答群の中から選べ。
（行番号）
01 function funcfdfx(x,n,a) 02 real funcfdfx 03 real,intent(in):: x 04 integer,intent(in):: n 05 real,dimension(0:10),intent(in):: a 06 real:: fx,dfx 07 integer:: k 08 fx = a(n) 09 dfx = real(n)*a(n) 10 do k = n-1,0,-1 11 fx = fx*x+a(k) 12 if(k>0) then 13 dfx = dfx*x+real(k)*a(k) 14 end if 15 end do 16 funcfdfx = fx/dfx 17 return 18 end function funcfdfx
解答群
ア fd = funcfdfx(xk,n,a)
イ fd = funcfdfx(xk,n-1,a)
ウ fd = xk-funcfdfx(xk,n,a)
エ fd = xk-funcfdfx(xk,n-1,a)
オ fd = xk-funcfdfx(xk1,n,a)
カ fd = xk-funcfdfx(xk1,n+1,a)
設問３ このプログラムと実行結果に関する次の記述中の に入れる正しい答えを，解答群の中から選べ。
このプログラムを用いて，方程式
f (x) ＝ x3 −7.300672 x2 ＋ 14.64975 x − 8.944963 ＝ 0
k x(k+1) x(k) f(x(k))/f'(x(k)) 1 1.1954747 1.0000000 -0.1954747 2 1.3013405 1.1954747 -0.1058657 3 1.3568242 1.3013405 -0.0554837 4 1.3852991 1.3568242 -0.0284749 5 1.3997463 1.3852991 -0.0144472 6 1.4070501 1.3997463 -0.0073039 7 1.4107746 1.4070501 -0.0037245 8 1.4127800 1.4107746 -0.0020054 9 1.4141418 1.4127800 -0.0013617 10 1.4215767 1.4141418 -0.0074349 no convergence
なお，この方程式の三つの解は，1.414213 …，1.414324 …，4.472135 … である。
この出力結果において，x(k+1) の値は，反復回数 k が 9 までは解に近づいているが，10 回目の計算で解から離れて収束していない。
収束しなかった原因を調べるために，24 行目と 25 行目の間に xk1，fx，dfx の値を出力する文を挿入して実行したところ，次の出力が得られた。
xk1 fx dfx 1.1954747 -0.5958862 3.0484056 1.3013405 -0.1568604 1.4816923 1.3568242 -0.0404415 0.7288904 1.3852991 -0.0102854 0.3612089 1.3997463 -0.0025959 0.1796818 1.4070501 -0.0006533 0.0894413 1.4107746 -0.0001650 0.0442972 1.4127800 -0.0000429 0.0213995 1.4141418 -0.0000124 0.0091047 1.4215767 -0.0000057 0.0007696
この結果が示すように，xk1 の値が解に近づくに従って，dfx の値は小さくなっている。したがって，収束しなかった理由として，二つの解の値が近いために，fx/dfxの値が ことによって，abs(xk1-xk) の値が ことがあげられる。
a に関する解答群
ア xk に近づかなかった イ xk に近づいた
ウ 0 に近づいた エ 0 に近づかなかった
b に関する解答群
ア eps より大きくならなかった イ eps より小さくならなかった
ウ eps と等しくなった エ 0 になった