質問をすることでしか得られない、回答やアドバイスがある。

15分調べてもわからないことは、質問しよう!

新規登録して質問してみよう
ただいま回答率
85.50%
FORTRAN

FORTRAN(フォートラン)は科学時術計算に向いた手続き型プログラミング言語です。 並列計算の最適化が行いやすい特性上、数値予報および気候モデルなどの大規模な計算を行う分野のスーパーコンピュータで使われています。

Q&A

解決済

1回答

3192閲覧

Fortran: べき級数の和を特定次とその前の項で比較

退会済みユーザー

退会済みユーザー

総合スコア0

FORTRAN

FORTRAN(フォートラン)は科学時術計算に向いた手続き型プログラミング言語です。 並列計算の最適化が行いやすい特性上、数値予報および気候モデルなどの大規模な計算を行う分野のスーパーコンピュータで使われています。

0グッド

0クリップ

投稿2017/03/13 14:02

編集2017/03/13 17:54

###前提・実現したいこと
現在,fortranの教科書を独学しています。問題を解きながら進めているのですが、どうしても組み方がわからない問題があるので教えてください。

問題は画像のようなものです。第n項目を出力することはできるのですが、n項目とn-1項目を比較しつつそれらの誤差を出力する方法がどうしてもわかりません。nとn項目とn-1項目と誤差(書籍での|e|)をファイルに出力したいです.

![問題内容:

以下に挙げたソースコードは、とりあえず一つの値を出力するものです。
(k=0.64ではK(0.64)=1.78423 63259 46738...)
####ソースコード(Fortran)
program first_eliptic_integral
implicit none
integer :: i, n=50 , fo=11 !50次まで足し上げ
real(8) :: pi=acos(-1.0d0), k, K1, X, sum=1.0d0, sum0 = 1.0d0

open(fo, file='Ex1.23_output.d') write(*,*) 'k(0<=k^2<=1):' read(*,*) k do i=1,n X=(k**2*(2*i-1)**2)/((2*i)**2) sum0 = sum0*X sum = sum + sum0 K1= (pi/2.0d0)*sum enddo write(fo,*) K1 close(fo)

end program first_eliptic_integral

###試したこと
nをカウントしたいので、do i=1,n文の外に更にnをカウントするdo文を置いてみたり、if文を使ってn-1を取り出そうとしましたが、出力結果は全部同じ数値になってしまいます。

###補足情報(言語/FW/ツール等のバージョンなど)
現在配列は勉強中なので、do文, if文(exit, cycle)のみで教えていただけたらうれしいです。(GOTO文もなしでお願いします.)

参考書籍は数値計算のためのFortran90/95プログラミング入門(牛島)です。

環境はgfortranでコンパイルしています。

※ここを利用するのは初めてなので、よろしくお願いします。

気になる質問をクリップする

クリップした質問は、後からいつでもMYページで確認できます。

またクリップした質問に回答があった際、通知やメールを受け取ることができます。

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

guest

回答1

0

ベストアンサー

今プログラム中の pi/2 * sum0 が第N項にあたるわけだから、それをループの内側に書けばいいのではないかと思う。

write(*, *) i, pi/2.0d0 * sum0

投稿2017/03/13 17:43

curehoney

総合スコア249

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

退会済みユーザー

退会済みユーザー

2017/03/13 17:52 編集

回答ありがとうございます。 ご指示を参考にしてdoループ内を以下のように書き換えてみました。 do i=1,n X=(k**2*(2*i-1)**2)/((2*i)**2) sum0 = sum0*X sum = sum + sum0 K2= (pi/2.0d0)*sum0 K1= (pi/2.0d0)*sum e=abs(K1-K2) write(fo,*) i, K1, K2, e enddo ところが結果は(i, K1, K2, e) 1 1.7316458706586941 0.16084954386379741 1.5707963267948968 2 1.7687056055649129 3.7059734906218930E-002 1.7316458706586939 3 1.7792470412715709 1.0541435706657828E-002 1.7687056055649131 4 1.7825528355091789 3.3057942376078947E-003 1.7792470412715711 5 1.7836496186981554 1.0967831889765968E-003 1.7825528355091789 6 1.7840271070988412 3.7748840068598958E-004 1.7836496186981552 7 1.7841604267573499 1.3331965850839719E-004 1.7840271070988414 8 1.7842084218344127 4.7995077063022986E-005 1.7841604267573496 ...以下50項目まで続く。 のようになり、eが0に収束しないようです。n-1項目をどこで取り出すかが難しいです。 もしご指示の内容を誤解している箇所がございましたらまた指摘していただけるとうれしいです。 ありがとうございます。
curehoney

2017/03/14 02:49

倍精度の有効桁数は14~15桁なので、1程度の量に足す項として1d-14程度以下になれば十分収束したと言えると思います。相対誤差として |足す項/総和| を取ってそれを見るのが一般的です。
退会済みユーザー

退会済みユーザー

2017/03/14 13:04 編集

コメントありがとうございます! ご指示いただいたとおり、単純な差ではなく相対誤差をとって以下のように書き換えてみました。 do i=1,n X=(k**2*(2*i-1)**2)/((2*i)**2) sum0 = sum0*X sum = sum + sum0 K2= (pi/2.0d0)*sum0 K1= (pi/2.0d0)*sum e=K2/K1 write(fo,*) i, K1, K2, e enddo この結果が 1 1.7316458706586941 0.16084954386379741 9.2888243831640044E-002 2 1.7687056055649129 3.7059734906218930E-002 2.0953026207197606E-002 3 1.7792470412715709 1.0541435706657828E-002 5.9246610853567593E-003 4 1.7825528355091789 3.3057942376078947E-003 1.8545280519910133E-003 5 1.7836496186981554 1.0967831889765968E-003 6.1490955257070801E-004 6 1.7840271070988412 3.7748840068598958E-004 2.1159342208642543E-004 7 1.7841604267573499 1.3331965850839719E-004 7.4724030703170053E-005 8 1.7842084218344127 4.7995077063022986E-005 2.6899927427580135E-005 9 1.7842259569839505 1.7535149537929349E-005 9.8278749220590346E-006 10 1.7842324390974693 6.4821135187891140E-006 3.6329983564630214E-006 となりまして、gnuplotで誤差を表示させても十分0に収束していると考えて良いことがわかりました。誤差の概念がまだ未熟な中、相対誤差というのは思いつきませんでした。ご教授いただきありがとうございます! 書籍とは少しやり方が違いますが、相対誤差を用いても同じ結果が導けるということを学びました。お付き合い頂きありがとうございました。
guest

あなたの回答

tips

太字

斜体

打ち消し線

見出し

引用テキストの挿入

コードの挿入

リンクの挿入

リストの挿入

番号リストの挿入

表の挿入

水平線の挿入

プレビュー

15分調べてもわからないことは
teratailで質問しよう!

ただいまの回答率
85.50%

質問をまとめることで
思考を整理して素早く解決

テンプレート機能で
簡単に質問をまとめる

質問する

関連した質問