くだすれFORTRAN（超初心者用）

1 名前：デフォルトの名無しさん [2006/01/24(火) 09:48:23 ] このスレッドは、他のスレッドでは書き込めない超低レベル、  もしくは質問者自身何が何だが分からない質問を勇気を持って書き込むスレッドです。  FORTRAN使いが優しくコメントを返しますが、  お礼はFORTRANの布教と初心者の救済をお願いします。  610 名前：デフォルトの名無しさん [2006/07/07(金) 11:11:35 ] fortranでC言語の#defineみたいに、プリプロセッサで定数指定できるのってありますか？ 611 名前：デフォルトの名無しさん mailto:sage [2006/07/07(金) 13:02:41 ] 拡張子を .F にしたらプリプロセッサを解析してくれるコンパイラも多い。 612 名前：610 [2006/07/07(金) 13:35:41 ] >>611 返信どうもっす。じゃあ、fortranの標準では#defineみたいなのないんですね？ 613 名前：デフォルトの名無しさん [2006/07/07(金) 13:56:35 ] >>611 標準では存在しない。 #define は、便利な反面、批判も多いのでたぶん意図的に新しい規格でも不採用とおもわれ。 ただ問題点を無くした条件コンパイルなどのプリプロセッサへの対応は論じつづけられている。 614 名前：610 [2006/07/07(金) 18:17:48 ] >>613 なるほど、納得いたしました。どうもありがとうございます！ 615 名前：デフォルトの名無しさん mailto:sage [2006/07/07(金) 20:37:06 ] 定数は parameter だな。 616 名前：デフォルトの名無しさん [2006/07/11(火) 19:31:21 ] 学校の課題をやっていただきたいのですが、、、 次の数列の和を求めよ。（DO文） 1/1+x=1-x+x**2-x**3+x**4-x**5・・・ 打ち切り条件　1.0*10**-6 提出明日です。よろしくおねがいします。 617 名前：デフォルトの名無しさん mailto:sage [2006/07/11(火) 22:20:20 ] 割り込みすいません… 分数の足し算プログラムですが、 どこが違うかわからないので教えてくださいませんか？ PROGRAM FRACTION INTEGER I,J,K,L,M,N READ(*,*) I,J,K,L CALL TUUBUN(I,J,K,L) WRITE(*,*)M,N CALL REDUCTION(M,N) M=M/A N=N/A WRITE(*,*)M,'/',N END SUBROUTINE TUUBUN(A,B,C,D) INTEGER A,B,C,D M=A*D+B*C N=C*D RETURN END 618 名前：デフォルトの名無しさん [2006/07/11(火) 22:24:36 ] >>616 FORTRAN７７で適当に作った。適当に直してくれ。 数学的には　-1<x<1 で収束する式だが、１に近づくにつれて収束が遅くなる。 つまり、たくさんの項の和を取らねばならない。しかし、延々とやると数値誤差が たまるので意味がなくなる。ゆえに適当なところであきらめて打ち切る。 PROGRAM VIP PARAMETER (eps = 1.0e-6) WRITE(6, *) 'INPUT x (-1 < x < 1)' READ(5, *) x s = 1.0 t = 1.0 DO 10 i = 1, 10000 t = - t * x s = s + t IF (ABS(t) .LT. eps) GOTO 999 10 CONTINUE STOP 'NOT CONVERGED' 999 WRITE(6, *) ' sum =', s, ' iteration', i WRITE(6, *) ' 1/(1+x) =', 1.0 / (1.0 + x) STOP END 619 名前：デフォルトの名無しさん [2006/07/11(火) 22:32:24 ] >>617 サブルーチンの変数の範囲について根本的に誤解している。 メインプログラムとサブルーチンでは変数はお互いに見えない。 ゆえに、通分した分子と分母はM、Nはサブルーチンの引数で返す必要がある。 CALL TUUBUN(I,J,K,L, m, n) SUBROUTINE TUUBUN(iA,iB,iC,iD, m, n) M=iA*iD+iB*iC N=iC*iD RETURN END あと、暗黙の型は死んでも守るべし。周りがなんと言おうと守るべし。 かっぱに尻子玉を抜かれそうになっても寧ろ暗黙の型守るべし。 620 名前：617 mailto:sage [2006/07/12(水) 10:25:30 ] >>619 途中からなくなっていたので、 もう一度直して PROGRAM FRACTION INTEGER I,J,K,L,M,N READ(*,*) I,J,K,L CALL TUUBUN(I,J,K,L) CALL REDUCTION(M,N) M=M/A N=N/A WRITE(*,*)M,'/',N END SUBROUTINE TUUBUN(iA,iB,iC,iD) INTEGER iA,iB,iC,iD M=iA*iD+iB*iC N=iC*iD RETURN END 621 名前：617 mailto:sage [2006/07/12(水) 10:26:19 ] >>620 SUBROUTINE REDUCTION(P,Q) INTEGER P,Q,X,Y X=P Y=Q 10 IF(X.EQ.Y) THEN GO TO 20 ELSE IF(X.GT.Y) THEN X=X-Y ELSE Y=Y-X END IF END IF GO TO 10 20 CONTINUE A=X RETURN END 622 名前：デフォルトの名無しさん [2006/07/12(水) 14:25:40 ] >>621 PROGRAM FRACTION READ(*,*) I, J, K, L CALL TUUBUN(I, J, K, L, M, N) CALL REDUCTION(M, N, IA) M = M / IA N = N / IA WRITE(*,*) M, '/', N END SUBROUTINE TUUBUN(N1, N2, ID1, ID2, NUMERA, IDENOM) NUMERA = N1 * ID2 + N2 * ID1 IDENOM = ID1 * ID2 RETURN END SUBROUTINE REDUCTION(NUMERA, IDENOM, IGCD) IX = NUMERA IY = IDENOM 10 IF(IX .EQ. IY) THEN GO TO 20 ELSE IF (IX .GT. IY) THEN IX = IX - IY ELSE IY = IY - IX END IF GO TO 10 20 CONTINUE IGCD = IX RETURN END >>622 623 名前：さおり [2006/07/12(水) 16:17:47 ] すいません！無理を承知でカキコします。 私は今理系で某国立大学に通ってるんですけど、水商売のバイトにハマってしまって全然授業についていってない状態なんです。 それでプログラミングの授業があるんですけど課題が出て「放物の運動（空気抵抗を考慮する）」のプログラムをfortranで作れって課題が出たんです。。 初速度と投げる角度を入力して一定時間ごと（１秒ごととか）の物体の位置をｘとｙで出すらしいんです。 私、授業にほとんど出てなくて本当にまったくわかってない状態なんです。 少し真面目な大学なんで私みたいな子は浮いてしまって、学校には助けてくれる友達もいないんです。 誰かプログラム作ってください。お願いします。この授業を落としてしまうと本当に留年なんです。おねがいします！！ 624 名前：デフォルトの名無しさん [2006/07/12(水) 16:43:42 ] >>623 まずはおっぱいうｐ汁 625 名前：デフォルトの名無しさん mailto:sage [2006/07/12(水) 21:23:23 ] >>623 担当教官をお店で接待すれば優もらえるおｗ 626 名前：デフォルトの名無しさん mailto:sage [2006/07/13(木) 00:14:44 ] 金があるなら留年しろ。 627 名前：デフォルトの名無しさん [2006/07/13(木) 01:04:04 ] 俺は少年なら助けてやるのだが・・・・女は興味ないｗｗｗ 628 名前：デフォルトの名無しさん mailto:sage [2006/07/13(木) 07:18:28 ] >>623 課題だけ書く方が答えが得られるすれのようです。 629 名前：さおり [2006/07/13(木) 10:47:04 ] 本当にお願いします！ いろいろ事情があって留年できないんです。 630 名前：デフォルトの名無しさん mailto:sage [2006/07/13(木) 12:22:06 ] >>629 課題は省略しないで出題されたまま提示するべき ヒントになるかどうか知らん ttp://www12.plala.or.jp/ksp/mechanics/resistdown/ 631 名前：デフォルトの名無しさん mailto:sage [2006/07/13(木) 12:26:25 ] >>623 空気定数って物体のサイズによって二通りの計算方法があるけど、どっちを使うの？ 632 名前：631 mailto:sage [2006/07/13(木) 12:29:35 ] >>629 >>630 のヒントの方の計算方法だと楽でいいよね。さあどっちだろうか？ 633 名前：631 mailto:sage [2006/07/13(木) 16:25:24 ] >>631 ×空気定数→○空気抵抗 634 名前：さおり [2006/07/13(木) 21:43:35 ] えーっと空気抵抗が速度の２乗に比例する場合らしいです。本当にありがとうございます。 635 名前：デフォルトの名無しさん mailto:sage [2006/07/13(木) 22:08:31 ] 解決したみたいだから、次の宿題。 636 名前：デフォルトの名無しさん [2006/07/13(木) 22:56:22 ] これから宿題をストーリーつきで頼む人は、 『ボク宿題が解けなくて困ってるんだ。 お兄ちゃん様に一緒に解いてほしいなぁ〜。云々』 という、もっとやる気を出させるストーリーを考えること！ 637 名前：デフォルトの名無しさん mailto:sage [2006/07/13(木) 23:00:37 ] >>634 ああ、めんどくさい方ね。 それ、解析的に解が求まらないから、微分方程式出して数値解析するしかないかな？ 一応、微分方程式を立ててみるけど、チェックしてみて vx, vy を速度の x,y 成分として、物体の質量を m 、重力定数を g、空気抵抗の 係数を k とすると m dvx/dt = -k vx √(vx^2+vy^2) m dvy/dt = -k vy √(vx^2+vy^2) - mg となるかな？あとは、授業で習った数値解析で vx, vy をまず求めるのだが…… 638 名前：デフォルトの名無しさん mailto:sage [2006/07/13(木) 23:19:43 ] >>636 兄貴ぃ〜！ 俺、宿題分かんないッス！！ 徹夜で手取り足取り教えて欲しいッス！！ 639 名前：デフォルトの名無しさん [2006/07/14(金) 00:15:18 ] >>638 サブ！　手取り腰取り教えてやるぜ！！ 640 名前：デフォルトの名無しさん [2006/07/14(金) 01:05:54 ] >>637 最後まで解けよｗ 641 名前：デフォルトの名無しさん mailto:sage [2006/07/14(金) 02:10:46 ] >>637 言いだしっぺのh(ry 642 名前：641 mailto:sage [2006/07/14(金) 02:11:58 ] 間違った 637でなくて>>640 643 名前：637 mailto:sage [2006/07/14(金) 10:13:37 ] >>640 そうじゃなくて、学校で習った数値解析法で解かないと 合格になる気がしないんですけど。XX法とか言ってもらえれば 続きも出来るのですが…… 644 名前：デフォルトの名無しさん mailto:sage [2006/07/14(金) 11:08:41 ] >>634の書き込みを見ると、「らしい」なんて書いているから、課題の正確な内容を 把握してなさそう。仮に折角>>637が解いたところで「やっぱり違った」とかで別解を 要求してきそうだ。 645 名前：デフォルトの名無しさん [2006/07/14(金) 13:55:36 ] C++スレ見たいに宿題依頼のテンプレが必要なんでね？ 見たところ、宿題はいまだF77中心のようだし。ｆ９０かｆ７７かも選択してもらわないと。 646 名前：さおり [2006/07/14(金) 16:17:43 ] そんなこと絶対しません。よろしくおねがいします。 647 名前：デフォルトの名無しさん mailto:sage [2006/07/14(金) 18:27:20 ] つーことで、さおりタンの仕事は 1. 微分方程式の確認 2. 授業でならった微分方程式の数値解析法の(せめて)名前の確認 だね。 648 名前：デフォルトの名無しさん mailto:sage [2006/07/14(金) 19:44:38 ] >>646 「そんなことしません」じゃなくて、レポートとして期待されているプログラムを 作成するためには情報が足りねーの。 課題を入手してきなさい。あと、おっぱいうｐな。 649 名前：デフォルトの名無しさん [2006/07/14(金) 20:48:58 ] >>637に従って解いたｗ 微分方程式は、一番原始的な方法で積分したｗ　 大学１，２生ならこれでも合格ラインを超えるだろう。３，４年生ならルンゲ・クッタを使うべきだろう。 原点（０，０）から初速ｖｘ、ｖｙで玉を放り投げる。地面より下に座標がきたら（ｙ＜０．０）終了する。 PROGRAM yamajun PARAMETER( g = 9.8, fric = 0.3 ) t = 0.0 dt = 1.0 / 8.0 amass = 1.0 x = 0.0 y = 0.0 vx = 10.0 vy = 10.0 ax = 0.0 ay = 0.0 1 CONTINUE WRITE(9, *) x, y, t !, vx, vy !, ax, ay IF ( y < 0.0 ) GOTO 999 t = t + dt x = x + dt * vx y = y + dt * vy ax = - fric * vx * SQRT(vx**2 + vy**2) / amass ay = -g - fric * vy * SQRT(vx**2 + vy**2) / amass vx = vx + dt * ax vy = vy + dt * ay GOTO 1 999 STOP END FORTRANは、うほっの世界であるべきなので、さおりは次回からは カヲルという少年に性転換すること。 650 名前：637 mailto:sage [2006/07/14(金) 21:17:58 ] >>649 thx. いやー。数値解析の教科書とか見て気が重かったのだが、ありがとう。 さおりタン満足してくれるといいね。 651 名前：649 mailto:sage [2006/07/14(金) 21:40:27 ] >>650 >>649はあまりに手抜きｗ　誰かもっと良いの作ってやれｗ ２元２階の微分方程式なので、本のルンゲ・クッタルーチンにそのまま放り込めないｗ ちょっと脳みそ使う必要がある。数独１回分くらいｗ　 今日暑いから考えらんねー　冷房ほしいよー 652 名前：デフォルトの名無しさん [2006/07/14(金) 21:43:43 ] 2,3,5.6,2.7,2,-,9,3 という「-」のような記号の混じった数列データが入ったファイルがあって、 それを下のように読み込みたいのですが、 real*8 a,b,c,d,e,f,g,h read(10,*) a,b,c,d,e,f,g,h これで「-」を読み込むときエラーが出て終了します。 「-」が出てきた場合には-1.0d-10などの任意の数字を代入しておく ような方法はないでしょうか？ よろしくお願いします。 653 名前：デフォルトの名無しさん [2006/07/14(金) 23:37:38 ] >>652　ほい。　F９０で書いたが、主要なところはF７７でもおｋ。　質問あればおｋ。 やり方は色々あるが、基本的にはまずデータを１行文字列として読み込んで、 カンマで区切って１こづつ数値として読み取るのが標準だろう。 ここでは読み取りエラーが起きてから例外処理的に”−”を処理しているが、 ”−”が多いならはじめから一度文字列で読み込んで処理するほうがすっきりする。 PROGRAM vip IMPLICIT NONE CHARACTER (LEN = 136) :: text REAL(KIND = 8) :: x(8) INTEGER :: i, ipos, io DO READ(9, *, IOSTAT = io) x IF (io == -1) EXIT IF (io > 0) THEN BACKSPACE(9) READ(9, '(A)') text DO i = 1, 8 ipos = INDEX(text, ',') - 1 IF (ipos == -1) ipos = LEN(text) ! 8th number IF (TEXT(ipos:ipos) == '-') THEN x(i) = -1.0d-10 ELSE READ(TEXT(1:ipos), *) x(i) END IF text = text(ipos + 2:) END DO END IF PRINT '(8g9.2)', x END DO STOP END PROGRAM vip 654 名前：６５３ [2006/07/14(金) 23:39:00 ] ■入力データ 2,3,5.6,2.7,2,1,9,3 2,3,5.6,2.7,2,-,9,3 2,3,5.6,2.7,2,2,9,3 ■出力 2.0 3.0 5.6 2.7 2.0 1.0 9.0 3.0 2.0 3.0 5.6 2.7 2.0 -0.10E-09 9.0 3.0 2.0 3.0 5.6 2.7 2.0 2.0 9.0 3.0 Press any key to continue 655 名前：デフォルトの名無しさん [2006/07/15(土) 03:04:03 ] 　　　　　 |　|/ノ二__‐──ｧ　　 ヽﾆニ二二二ヾ } ,'⌒ヽ 　　　　 /⌒!|　 =彳ｏ｡ﾄ￣ヽ　　　　　'´ !o_ｼ｀ヾ　| i/ ヽ !　　　　ウホッ！　いい男・・・ 　　　　 ! ハ!|　　ー─　'　　i　　!　 　 `' 　 '' " 　 ||ヽ　l | 　　　　| | /ヽ!　　　　　　　　|　　　　　　　　　　　 |ヽ i　! 　 656 名前：652 mailto:sage [2006/07/15(土) 09:43:09 ] >>653 ウホッ！ 657 名前：デフォルトの名無しさん mailto:sage [2006/07/15(土) 09:44:53 ] >>623,>>629,>>634,>>646 >>637をもとに>>649が解いたものは、m、k、初期位置、時間の刻みに 適当な値を代入して解いている。本来はあんたが提示すべき情報。 >>630,>>647,>>648を良く読め。 658 名前：デフォルトの名無しさん mailto:sage [2006/07/15(土) 10:57:31 ] 教えて！goo にも同じような質問が・・・ 659 名前：637 mailto:sage [2006/07/15(土) 12:20:32 ] >>658 じゃあ我々の役目も終ったのかな？ みんなお疲れ様。 660 名前：デフォルトの名無しさん [2006/07/15(土) 17:35:26 ] 女が来るとふいんきが乱れるから、ホモOnlyスレにしないか？ 661 名前：デフォルトの名無しさん mailto:sage [2006/07/15(土) 18:59:43 ] プログラムの実行中に一瞬ポーズする関数ってありますか？ 662 名前：デフォルトの名無しさん mailto:sage [2006/07/15(土) 19:17:53 ] >>661 PAUSE 663 名前：デフォルトの名無しさん mailto:sage [2006/07/15(土) 19:30:41 ] >>662 ただPAUSEはFORTRAN2003あたりで廃止勧告が出されている。 TSS利用で便利だから漏れは好きなんだけどなｗ 664 名前：デフォルトの名無しさん [2006/07/15(土) 20:28:29 ] write(*,*) abs(-16.8-(6.8))>10.0 を出力すると「T」になってしまいます。 助けてください。 665 名前：664 [2006/07/15(土) 20:32:23 ] 間違えました、 real*8 a,b a=-16.80000000000000 b=- 6.80000000000000 write(*,*) abs(a-b)>10 とやったら「T」と出てしまいます。 何が原因でしょうか？ 666 名前：664 mailto:sage [2006/07/15(土) 20:53:41 ] 更に状況を書かせて頂きますと、 読み取りファイルに -16.8,-6.8 とあって、 real*8 a,b open(10,file="test") read(10,*)a,b write(*,*) abs(a-b)>10 とやったら「T」と出力されるのです。 667 名前：デフォルトの名無しさん mailto:sage [2006/07/15(土) 20:59:25 ] FORTRAN77は大文字でプログラム書くのがデフォルトなのか 668 名前：デフォルトの名無しさん [2006/07/15(土) 21:04:15 ] >>664 それは良くあることｗ　気にするな。　浮動小数点之罪。 669 名前：デフォルトの名無しさん [2006/07/15(土) 21:06:04 ] >>667 本来のFORTRAN77の規格では大文字のみ。 そもそも、昔は小文字が表示できない端末とか良くあったｗ 670 名前：デフォルトの名無しさん [2006/07/15(土) 21:22:32 ] >>664 参考までに、INTEL　FORTRANの結果を見てくれ。 PROGRAM unko IMPLICIT NONE REAL (KIND = 8) :: a, b a = -16.8 b = -6.8 PRINT *, a, b, ABS(a - b) < 10 a = -16.8d0 b = -6.8d0 PRINT *, a, b, ABS(a - b) < 10.0d0 STOP END PROGRAM unko ■結果 -16.7999992370605 -6.80000019073486 T -16.8000000000000 -6.80000000000000 F Press any key to continue この場合は倍精度（８バイト）の変数に単精度（４バイト）の常数を与えているので、 長さの足りない４バイトにゴミが入っている。 ■ とにかく浮動小数点の誤差とか桁落ちとか実数の比較とかでググレば説明が出るはず。 この辺は昔は電子計算機利用の一番最初に注意されて来た事柄だ。 671 名前：664 mailto:sage [2006/07/15(土) 21:45:51 ] >>670 どうもありがとうございます。 一応計算中の変数を画面に表示させてチェックしたのですが、 aは-16.8000000000000 bは -6.8000000000000 と表示されていました。 この上でabs(a-b)>10をやるとTになりました。 672 名前：デフォルトの名無しさん mailto:sage [2006/07/15(土) 22:03:21 ] >>671 (T_T) 30桁くらい出力してごらん。 write ( *, '(2f35.30)' ) a, b 673 名前：670 [2006/07/16(日) 00:28:13 ] >>671 漏れの書き方が誤解を招いたらしい。すまんこ。 数直線のうちFORTRANであらわされるのは、０を中心として正負にほぼ同じだけの 範囲の中の、有限個の有理数だけである。数直線の上に、定規の目のように切った点が あって、定規の目に挟まれた部分は近くの目で代表することになる。 （浮動小数点のばあい、この定規の目は体育で使う巻尺のように、０の近傍は細かく 区切ってあるが、遠くへ行くほど荒く区切るようになっている。） コンパイラにも依存するのだが、基本的には単精度や倍精度の浮動小数点は、 ２の（負を含む）べき乗の和であらわされるような有理数になっている。 このため、１０進表記で書いた有限小数が、かならずしも数直線の上に切った目に 一致するとは限らない。 つまり １．５＝１＋０．５＝２＾０＋２＾−１ ０．７５＝０．５＋０．２５＝２＾−１＋２＾−２ これらの数は、きっかり割り切れる（すなわち、定規の目の上にのっている） ところが ０．８＝０．５＋０．２５＋０．０６２５＋・・・・＝２＾−１＋２＾−２＋２＾−４＋・・・・・ となって、数直線の上に切った目の上に載っていない。 つまり、近くの別の数に四捨五入されている。 >>672のように、FORMATで小数点いかの数を長く出すと、以下のようになる。 -16.7999992370605 -6.80000019073486 T -16.8000000000000 -6.80000000000000 F -16.800000000000000710542735760100 -6.799999999999999822364316059975 つまり、倍精度の−１６．８は、実際は-16.800000000000000710542735760100・・・・・・ という数に四捨五入されている。−６．８も同様。 （ただし、有効桁からすれば、＊フォーマットで出力される分しか制度が無い。 つまり＊フォーマットで出力される数の最後の桁程度の範囲で四捨五入がされている。） 674 名前：デフォルトの名無しさん [2006/07/16(日) 14:38:42 ] 　積分区間の下限Aと上限B、分割数N（偶数）を読み込んで、 Simpson公式を用いて、定積分I=(A,B)∫f(x)dxを計算するプログラムを作れ。 　ただし、数値積分の部分のプログラムは、任意の関数に対して使用できる サブルーチン・サブプログラムを用いよ。その場合関数f(x)の計算は関数サブプログラムで行う。 　さらに上記の分割数Nが奇数に対しても使えるように拡張せよ。 Simpson公式を用いて、Nが奇数のときは1区画だけ台形公式を用いるようにすればよい。 手も足も出ません…、お願いできますか？ 675 名前：デフォルトの名無しさん mailto:sage [2006/07/16(日) 16:37:20 ] >>674 PROGRAM HOGE READ(*,*) A, B, N IF (MOD(N,2) .NE. 0) GOTO 9999　… ※１ C DX = (B - A) / N S = 0. DO 10 I = 1, N/2 X1 = A + 2*(I-1)*DX X2 = A + (2*I - 1)*DX X3 = A + 2*I*DX S = S + (F(X1) + 4*F(X2) + F(X3)) 10 CONTINUE S = S * DX / 3. C IF (MOD(N,2) .NE. 0) THEN　… ※２ C X1 = B - DX C X2 = B C S = S + (F(X1) + F(X2)) * DX / 2 C END IF WRITE(*,*) S 9999 CONTINUE STOP END 前半部分はそのまま。（※２の所はいらない） 後半部分は※１のIF文を消して※２のコメントアウトされてるIFブロックを追加する。 676 名前：デフォルトの名無しさん [2006/07/16(日) 21:55:10 ] ファイルを読み込んだり作ったりするとき、 実行ファイルと同じディレクトリの中に 別のディレクトリを作って保存したりするにはどうしたら良いでしょうか。 677 名前：デフォルトの名無しさん [2006/07/16(日) 22:07:35 ] >>676 既出の質問だが・・・・お答えしよう。 FORTRANの標準命令ではできない。 コンパイラメーカー固有の拡張命令でシステムの関数を呼び出すものがあるだろうから それらを通じてやるしかない。マニュアルの後ろのほうの章を探して味噌。 678 名前：676 mailto:sage [2006/07/16(日) 23:32:02 ] >>677 ｻﾝｸｽｺ 679 名前：デフォルトの名無しさん mailto:sage [2006/07/17(月) 12:21:13 ] 番組の途中ですが・・・・・・・・・・・・ 盧武鉉（ノ・ムヒョン）大統領が今月11日に行われた与党ヨルリン・ウリ党の指導部 および国会の統一外交通商委員会に所属する議員らとの晩さん懇談会で行った 発言が波紋を呼んでいる。 一部新聞は懇談会出席者の証言を引用し、盧大統領は「ブッシュ米大統領が 北朝鮮問題を善と悪の対立概念で見ているため、説得が難しくなっている。 米国は友邦なので厳しく責めることは出来ないが、 日本とは対決しなければならない」と語ったという。 ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝ ソース:朝鮮日報 japanese.chosun.com/site/data/html_dir/2006/07/17/20060717000017.html 680 名前：デフォルトの名無しさん [2006/07/18(火) 05:15:28 ] 多変数の4次のルンゲクッタのプログラムを書いているのですが、どうもうまくいきません。 今、厳密解を y(1) = exp(x) + exp(-x) y(2) = exp(x) - exp(-x) とおいて、その微分方程式 dydx(1) = exp(x) - exp(-x) dydx(2) = exp(x) + exp(-x) を数値的に解いて、厳密解と比較しています。 サブルーチンのderive中で、 dydx(1) = exp(x) - exp(-x) dydx(2) = exp(x) + exp(-x) とおいた場合は、きちんと数値解と厳密解が一致するのですが、 dydx(1) = y(2) dydx(2) = y(1) とおくと、両者は一致しません。これはなぜでしょうか？ 考えているのですが、一向にわかりません。 プログラムは以下になります。よろしくお願いします。 681 名前：デフォルトの名無しさん [2006/07/18(火) 05:17:34 ] program rk4 implicit none double precision x0, y0, x, h, y(2), dydx(2), yout(2) external derive h = 1.0d-2 x0 = 0.0d0 y(1) = 2.0d0 y(2) = 0.0d0 do x = x0, 2.0d0, h write(*, *) x, dabs((dexp(x) + dexp(-x)) - y(1)) write(*, *) x, dabs((dexp(x) - dexp(-x)) - y(2)) write(*, *) call derive(x, y, dydx) call rk4(y, dydx, x, h, yout, derive) y(1) = yout(1) y(2) = yout(2) end do stop end subroutine derive(x, y, dydx) implicit none double precision x, y(2), dydx(2) c dydx(1) = y(2) c dydx(2) = y(1) dydx(1) = dexp(x) - dexp(-x) dydx(2) = dexp(x) + dexp(-x) return end 682 名前：デフォルトの名無しさん [2006/07/18(火) 05:18:07 ] subroutine rk4(y, dydx, x, h, yout, derive) implicit none integer i double precision h, x, dydx(2), y(2), yout(2) external derive double precision h6, hh, xh, dym(2), dyt(2), yt(2) hh = h*0.5d0 h6 = h/6d0 xh = x + hh do i = 1, 2 yt(i) = y(i) + hh*dydx(i) call derive(xh, yt, dyt) yt(i) = y(i) + hh*dyt(i) call derive(xh, yt, dym) yt(i) = y(i) + h*dym(i) dym(i) = dyt(i) + dym(i) call derive(x+h, yt, dyt) yout(i) = y(i) + h6*(dydx(i) + dyt(i) + 2.0d0*dym(i)) end do return end 683 名前：デフォルトの名無しさん [2006/07/18(火) 15:03:58 ] >>680 うまく行かないのは、subroutine rk4 に誤りがあるから。 以下のループ部分は、i　=　1 について積分し、次にi　=　2について積分しているが、 これはおかしい。i=1,2は独立な量ではなく、相互に依存している。 i=1,2を同時に積分して、ちょっとづつ前進させなければならない。 do i = 1, 2 yt(i) = y(i) + hh*dydx(i) call derive(xh, yt, dyt) yt(i) = y(i) + hh*dyt(i) call derive(xh, yt, dym) yt(i) = y(i) + h*dym(i) dym(i) = dyt(i) + dym(i) call derive(x+h, yt, dyt) yout(i) = y(i) + h6*(dydx(i) + dyt(i) + 2.0d0*dym(i)) end do こう直せばおｋ yt = y + hh*dydx call derive(xh, yt, dyt) yt = y + hh*dyt call derive(xh, yt, dym) yt = y + h*dym dym = dyt + dym call derive(x+h, yt, dyt) yout = y + h6*(dydx + dyt + 2.0d0*dym) 数値誤差は１０＾−１０程度になるが、これは数値微分を使っているのだからまぁおｋ 684 名前：デフォルトの名無しさん [2006/07/18(火) 21:45:16 ] 683さん、お答え頂きありがとうございます！ おっしゃるとおり、サブルーチンを変えたところうまく回りました！ どうもありがとうございました！ いまいち、サブルーチンと配列の組み合わせが理解できないでいます。 例えば、yを配列として(例えばy(3)などと)定義した場合、 yt = y + hh*dydx で全てのyを計算するのと、 do i = 1, 2 yt(i) = y(i) + hh*dydx(i) で、ひとつずつyを計算するのは違うのでしょうか？ というのも、ルンゲクッタギルというサブルーチンをひろってきたのですが、 そのコードでは、iを1からnまで計算していています。これを参考に、 間違いを指摘して頂いたコードを作ったのですが・・・。 一体何が違うのでしょうか？お手を煩わせてすみません。 よろしくお願いします。 685 名前：デフォルトの名無しさん [2006/07/18(火) 21:46:09 ] subroutine rkg(derive, x0, y, n, ndim, h, nx, f) implicit real*8 (a-h, o-z) external derive dimension y(ndim, 0:*), f(ndim, 4) data fp / 1. 70710 67811 86547 52440 d0 / data fm / 0. 29289 32188 13452 47560 d0 / data fn / 0. 16666 66666 66666 66667 d0 / ! 1/6 data fp2 / 3. 41421 35623 73095 04880 d0 / ! 2 + sqrt(2) data fm2 / 0. 58578 64376 26904 95119 d0 / data fpp / -4. 12132 03435 59642 57320 d0 / data fmm / 0. 12132 03435 59642 57320 d0 / hh = 0.5d0*h fmh = fm*h fph = fp*h fnh = fn*h 686 名前：デフォルトの名無しさん [2006/07/18(火) 21:47:18 ] do 100 j = 1, nx x = dble(j-1)*h + x0 call derivs(x, y(1, j-1), f(1, 3)) do i = 1, n f(i, 2) = y(i, j-1) + hh*f(i, 3) end do call derivs(x + hh, f(1, 2), f(1, 4)) do i = 1, n f(i, 1) = f(i, 2) + fmh*(f(i, 4)-f(i, 3)) f(i, 2) = fm2*f(i, 4) + fmm*f(i, 3) end do call derivs(x + hh, f(1, 1), f(1, 3)) do i = 1, n f(i, 4) = f(i, 1) + fph*(f(i, 3)-f(i, 2)) f(i, 1) = fp2*f(i, 3) + fpp*f(i, 2) end do call derivs(x+h, f(1, 4), f(1, 2)) do i = 1, n y(i, j) = f(i, 4) + fnh*(f(i, 2)-2.0d0*f(i, 1)) end do 100 continue return end 687 名前：デフォルトの名無しさん [2006/07/18(火) 22:58:38 ] >>684 DO　LOOPを消したのは、めんどくさかったからで、あらわに書けば>>685のようになる。 >>685のコードと、元の>>682のコードを見比べてみるとわかるが、 >>685ではｄｅｒｉｖｅを前ごとにループを廻しているのに対して、>>682では全体に対して ループをまわしている。これでは意味が違ってしまう。 そもそも、もっとも素朴なオイラー法ではy(x+h)=y(x)+hy'(x)なのだが、誤差が大きい。 そこでルンゲ・クッタでは、一気にｈだけ進まないで、まずh/2進んで、そこでのｙをもとめて、 これを元にｙ’をもとめて、さらにまたh/2進むと言う感じになっている。 イメージ的には、オイラーは両足そろえてｈだけジャンプするが、ルンゲクッタは 右足１歩左足一歩、もう一回右足一歩、左足一歩だ。 ところが、>>682のように一気にDOLOOPをまわすと右足ケンケン２回、左足ケンケン２回 的な状況になっていて、なんだかおかしいのだ。 なっている。 （あくまでイメージなので、真面目に考えてくださいｗ） 688 名前：687 [2006/07/18(火) 23:06:43 ] 日本語狂いまくりんぐｗ do i = 1, 2 　yt(i) = y(i) + hh*dydx(i) end do call derive(xh, yt, dyt) do i = 1, 2 　yt(i) = y(i) + hh*dyt(i) end do call derive(xh, yt, dym) do i = 1, 2 　yt(i) = y(i) + h*dym(i) end do do i = 1, 2 　dym(i) = dyt(i) + dym(i) end do call derive(x+h, yt, dyt) do i = 1, 2 　yout(i) = y(i) + h6*(dydx(i) + dyt(i) + 2.0d0*dym(i)) end do こういう風に、個々にDoLoopになってればおｋ！　F90なら、DoLoop省略できる。 689 名前：デフォルトの名無しさん [2006/07/19(水) 07:56:54 ] こんにちは、fortran90での課題が手に負えないので助けてください。 「n×nの実数正方行列Ａの行列式｜Ａ｜を再帰手続きでもとめよ」というもの ですが、再帰手続きの使い方、使いどころが見当もつきません。３×３までは できたんですが再帰手続きを使っておらず、それ以上に拡張できません。 長いこと詰まっています。どうか教えてください。 690 名前：デフォルトの名無しさん mailto:sage [2006/07/19(水) 09:40:39 ] 再帰手続きでってことなら、 一番泥臭い方法を使えって事かな。 先ず、行列式は 　　　　　　　　　 p　　n ｜Ａ｜ = Σ(-1)　P Πa_ii 　　　　　　P　　　　　 i と書ける。 多分、線形代数の教科書見たら書いてあると思う。 ここで、Π_i^n a_ii は a_11 a_22 ... a_nn という積を表している。 で、P は添字の置換演算子。 a_11 a_22 ... a_nn を a_21 a_12 ... a_nn とか a_11 a_n2 ... a_2n とか、 行列の添字の一方に注目して、それを置換する演算子。 Σ_P は、あらゆる置換を考慮して、全部足せ、ということを表している。 で、p は偶置換だと 0 で、奇置換だと 1 になる。 この「あらゆる置換を考慮して」というところを、 再帰を使って実現すればいい。 691 名前：デフォルトの名無しさん [2006/07/19(水) 10:48:05 ] 690さん、ありがとうございます。アドバイスをもとに考えてみます。 fortranをはじめたばかりの初心者なので、不躾なんですがよろしければ なにか具体的な例を示していただければ幸いです。行列を読み込んで添字に 着目して、部分行列を作ってそれから再帰手続きとしようとしているんですが、 再帰手続きの手順がこんがらがりよく理解できません。どうかお願いします。 692 名前：デフォルトの名無しさん mailto:sage [2006/07/19(水) 11:19:36 ] >>689 690 ではない者ですが、余因子展開すると n×n 行列の行列式を n-1×n-1 行列の行列式で表せると思います。 693 名前：デフォルトの名無しさん [2006/07/19(水) 15:39:04 ] >>689 MODULE mod_det IMPLICIT NONE CONTAINS RECURSIVE FUNCTION det(a) RESULT(res) IMPLICIT NONE REAL (KIND = 8), INTENT(IN) :: a(:, :) REAL (KIND = 8) :: res, wk(SIZE(a, 1) - 1, SIZE(a, 1) - 1) INTEGER :: i, n n = SIZE(a, 1) res = 0.0d0 IF (n == 1) THEN res = a(1, 1) ELSE DO i = 1, n wk(:, 1:i - 1) = a(2:, 1:i - 1) wk(:, i:n - 1) = a(2:, i + 1: n) res = res + REAL((-1)**(i - 1), KIND = 8) * a(1, i) * det(wk) END DO END IF RETURN END FUNCTION det END MODULE mod_det >>693 694 名前：デフォルトの名無しさん mailto:sage [2006/07/19(水) 15:40:57 ] PROGRAM vip USE mod_det IMPLICIT NONE INTEGER, PARAMETER :: n = 7 INTEGER :: i, j REAL (KIND = 8) :: chk, d, a(n, n) DO i = 1, n DO j = 1, n a(i, j) = 1.0d0 / REAL(i + j, KIND = 8) END DO END DO d = det(a) chk = 1.0d0 DO i = 1, n DO j = 1, n IF ( i > j ) chk = chk * REAL(i - j, KIND = 8) * REAL(i - j, KIND = 8) chk = chk / REAL( i + j, KIND = 8 ) END DO END DO PRINT *, d, d - chk STOP END PROGRAM vip >>694 余因子展開で解いた。　上の１行を元に展開している。 メインルーチンでは、Ｃａｕｃｈｙの定理でチェックしている。多分大丈夫？ 695 名前：デフォルトの名無しさん [2006/07/19(水) 15:45:01 ] ＞＞６９３の式では、小行列式を作るのに、あらわにワーク行列を作っているが、 ＦＯＲＴＲＡＮ９０なら、マスキングとかを作って、ワークをあらわに経由しないで、 小行列を引数として再帰関数を呼び出せるような気がする。 手元に紙マニュアルがないので、調べる気がしない。だれか頼む。 でもどっちにしろmask行列を作らないと駄目かな？ 意味ないかもｗ 696 名前：デフォルトの名無しさん [2006/07/19(水) 15:51:12 ] 687さん、早速お答え頂きどうもありがとうございます！ 全体と個々のstepとでloopを回すのでは全く違いますね。 ありがとうございました。 697 名前：デフォルトの名無しさん [2006/07/19(水) 16:24:45 ] 692さん、693さん、694さん、695さん、具体例およびアドバイスありがとうございます。 本当に助かります。早速これらのプログラムを参考に自分なりに考えてみます。 698 名前：デフォルトの名無しさん [2006/07/20(木) 00:14:15 ] すいません、fortran77で、 「sin(x)をマクローリン展開し、第ｎ項まで計算したものを、 ｘを0から180度まで30度ずつ変えて表示するプログラムを作成せよ。 但しsin(x)の計算はサブルーチンを使うこと」 という課題が出されたのですが、うまく作れません。 どうかヒントだけでも教えてもらえないでしょうか。 699 名前：デフォルトの名無しさん mailto:sage [2006/07/20(木) 03:03:56 ] >>698 多分マクローリン展開が間違ってるんじゃないかな？ 700 名前：デフォルトの名無しさん [2006/07/20(木) 03:24:06 ] >>698 PROGRAM oppai pi = 4.0 * ATAN(1.0) DO 10 i = 0, 180, 30 rad = pi * REAL(i) /180.0 CALL sine(10, rad, s) WRITE(6, *) ' SIN(', i, ') =', s, SIN(rad) 10 CONTINUE STOP END C SUBROUTINE sine(n, rad, s) s = 0.0 x = rad DO 10 i = 1, n s = s + x x = - x * rad * rad / REAL( 2 * i * (2 * i + 1) ) 10 CONTINUE RETURN END 701 名前：700 [2006/07/20(木) 11:34:13 ] 補足すると、SINをテイラー展開で求めるサブルーチンは、 １番目の引数が第何項まで足すかと言うもの、２番めはSIN（Θ）のΘ、３番目が結果の値。 貼り付けたソースでは１０項まで足しているが、単精度なのでこんなに項を足しても意味無い。 有効桁からより適切な項数を編み出してくれ。 SIN(x)=x-x^3/3!+x^5/5!-x^7/7!+........(-1)^((n - 1)/2)x^n/n! (n odd) 702 名前：デフォルトの名無しさん mailto:sage [2006/07/20(木) 12:58:47 ] >>700 SIN(x)=x-x^3/3!+x^5/5!-x^7/7! これは偶数次の展開が0になって消えているから7次(または8次)の項まで 展開したという理解で良かったっけ? 703 名前：700 [2006/07/20(木) 13:50:19 ] >>702 ういうい。 sin は奇関数だから、偶べきの項は出てこない。ゆえに、偶数次の項は消える。 704 名前：困ってます [2006/07/20(木) 16:08:00 ] 2次元のﾎﾟｱｿﾝ方程式 方程式 d2f/dx2 + d2f/dy2 = c 境界値　　f = 0 範囲 　 x= -1 〜 1 , y= 0 〜 1 定数 c= -1 (x<0) , c = 0 (x=0) , c = 1 (x>0) 離散式 -4*f(i,j)+f(i-1,j)+f(i+1,j)+f(i,j-1)+f(i,j+1)= dx*dx*c これをvisual Fortranを使ってやるのですが全くわかりません。。。 誰か出来る人いませんか？？？？ あとDOSで開いた後立体的に見たいんですけど よかったら貼り付けてください！！！ 705 名前：デフォルトの名無しさん [2006/07/20(木) 16:15:40 ] >>704 こんなの平均値を求めるのと変わらんぞｗ 本質は１行だｗ あっちのスレで解いてやったが、今後はこっちで継続せよ。 706 名前：困ってます [2006/07/20(木) 16:50:26 ] あれでやってみたんですが全くできません。。。 なんとか馬鹿でも分かるようにおしえてください！！！！ 707 名前：デフォルトの名無しさん [2006/07/20(木) 16:57:12 ] >>706 配列ｆの中に、答えが入っている。 これを出力してグラフを書けばいい。 　＿ /　　\ 　　　　\＿/ ⇒ｘ軸 　\＿/ ⇒ｙ軸　　ｘ＜０ 　＿ /　　\ ⇒ｙ軸　　ｘ＞０ こんな感じになっている。 708 名前：困ってます [2006/07/20(木) 17:10:13 ] ありがとうございます！！！ でもそれでもいまいち分かってない自分がいます。。。 プログラムをそのままコンパイル出来る形で教えていただけないでしょうか？？？？ 709 名前：困ってます [2006/07/20(木) 17:13:06 ] 出力したいんですがやり方がやっぱわかないです！！！！！ 710 名前：デフォルトの名無しさん [2006/07/20(木) 17:56:17 ] >>709 ごめ、元の奴括弧の範囲を継続行に直したところで、ずらしてしまっていた。 結果は大して変わらんが。訂正しておくｗ PROGRAM unko IMPLICIT NONE INTEGER, PARAMETER :: nx = 19, ny = 9 REAL :: f(-nx - 1:nx + 1, 0:ny + 1), ff(-nx:nx, ny), c(-nx:nx, ny), dx, dy INTEGER :: i, ix, iy dx = 1.0 / REAL(nx + 1) dy = 1.0 / REAL(ny + 1) c(-nx:-1, :) = -1.0 c( 0, :) = 0.0 c( 1:nx , :) = 1.0 f = 0.0 f(-nx - 1, :) = 0.0 f( nx + 1, :) = 0.0 f(:, 0) = 0.0 f(:, ny + 1) = 0.0 DO i = 1, 100 DO ix = -nx, nx DO iy = 1, ny ff(ix, iy) = ( -c(ix, iy) * dx * dy & + f(ix - 1, iy - 1) + f(ix - 1, iy + 1) + f(ix + 1, iy - 1) + f(ix + 1, iy + 1) ) / 4.0 END DO END DO f(-nx:nx, 1:ny) = ff END DO DO ix = -nx - 1, nx + 1 WRITE(9, '(1000es11.3)') dx * ix , f(ix, :) END DO STOP END PROGRAM unko

---

---

[ 続きを読む ] / [ 携帯版 ]

read.cgi ver5.27 [feat.BBS2 +1.6] / e.0.2 (02/09/03) / eucaly.net products.
担当:undef