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

1 名前：デフォルトの名無しさん [2006/01/24(火) 09:48:23 ] このスレッドは、他のスレッドでは書き込めない超低レベル、  もしくは質問者自身何が何だが分からない質問を勇気を持って書き込むスレッドです。  FORTRAN使いが優しくコメントを返しますが、  お礼はFORTRANの布教と初心者の救済をお願いします。  101 名前：９５ mailto:sage [2006/02/04(土) 19:40:40 ] オッス！おれ極右！ >>98氏解説�d 最近寒いので湯たんぽ代わりにノートパソコンを腹の上において寝ているお。 昨夜そういう寝そべった状態で布団の中で作ったのが>>95だお。 オラはVipperだから間違っていても謝罪も賠償もしないお！ 最小二乗法なんて普段使わないから、きょうの昼おっきしてからNumericalRecipesで 式を確かめといたおｗ　方式はちょっと違うが式おｋｗｗｗ　（15.2.4~6)あたり。 ttp://www.library.cornell.edu/nr/bookfpdf/f15-2.pdf FORTRAN77で書くと３倍くらいの長さになるとおもうお。 MATMULやTRANSPOSEはF90にしかない行列積と転置だお。 LOGを取ったのは線形化するためだお。 普通、最小二乗法は線形の式にしか使わないと思ったから、そうしてみたおｗ NumericalRecipesは、旧版のC言語版なら訳本が出ているお。 ttp://www.amazon.co.jp/exec/obidos/ASIN/4874085601/qid%3D1139049291/249-3009705-2437946 本のコードにはバグがちりばめられているので、本を打ち込むのはお勧めしないお。 金出して買うコードのほうはこまごまと修正されているお。 （＾ω＾）このくらいソースコードが短いと勉強する気にもなるだろ？ 102 名前：デフォルトの名無しさん mailto:sage [2006/02/05(日) 07:53:27 ] >>87 フォートランはわかんないからベーシック あってるかどうかもわかんない こんにちは 10 gosub 60 20 nannichi0=nannichi 30 gosub 60 40 print "日数は";nannichi-nannichi0 50 end 60 input "年,月,日",year,month,day 70 if month<3 then year=year-1:month=month+12 80 nannichi=year+int(year/4)-int(year/100)+int(year/400)+int((153*month-457)/5)+day 90 return 103 名前：デフォルトの名無しさん mailto:sage [2006/02/05(日) 11:38:27 ] >>102 解説頼む 104 名前：デフォルトの名無しさん mailto:sage [2006/02/05(日) 20:01:35 ] >>103 まちがいた。ちょと訂正。 80 nannichi=365*year+int(year/4)-int(year/100)+int(year/400)+int(306*(month+1)/10)+day-428 「これは Fairfield公式 といって、西暦1年1月1日からの経過日数を算出する式らしい。 　ただしグレゴリオ暦が採用されたのは1582年10月15日からだが。 　要するにこれを使って2つの日付の差を求めればよい。」 っておばあちゃんがいってた。 105 名前：94 mailto:sage [2006/02/05(日) 20:12:40 ] すいませんがこないだのプログラムを ＤＯ文で作っていただけませんか。 何回もお願いして、すいません。 106 名前：デフォルトの名無しさん mailto:sage [2006/02/05(日) 20:33:52 ] コンピューターおばあちゃん♪ コンピューターおばあちゃん♪ イエーイ　イエーイ　僕は大好きさ〜♪ 107 名前：デフォルトの名無しさん mailto:sage [2006/02/05(日) 21:07:36 ] >>104 把握したｗ ２月の２８日の例外を２月をお尻にすることで回避。 大の月と小の月の入れ替えを１ヶ月の平均日数を30.6日にすることで実現。 かなり絶妙な公式だ。ｵﾓｼﾛｽｗ だが最近流行の明示的なプログラミングから離れているがｗｗｗｗｗ 108 名前：デフォルトの名無しさん mailto:sage [2006/02/05(日) 21:54:46 ] >>105　世話が焼けるッス！　FORT７７で、素朴に書いといたからｗ PROGRAM vipper 　　　PARAMETER (n = 5) REAL x(n), y(n) 　　　DATA x /0.1 , 0.2 , 0.3 , 0.4 , 0.5 / DATA y /1.228, 1.005, 0.823, 0.674, 0.552/ 　　　sumx = 0.0 　　　sumy = 0.0 　　　sumxx = 0.0 　　　sumxy = 0.0 sum = REAL(n) DO 10 i = 1, n 　　　 sumx = sumx + x(i) sumy = sumy + LOG(y(i)) sumxx = sumxx + x(i) * x(i) sumxy = sumxy + x(i) * LOG(y(i)) 10 CONTINUE det = sum * sumxx - sumx * sumx alog = ( sumxx * sumy - sumx * sumxy ) / det b = ( sum * sumxy - sumx * sumy ) / det WRITE(6, *) ' a = ', EXP(alog), ' b = ', b 　　　END >>108　インデント 109 名前：やり直しｗ mailto:sage [2006/02/05(日) 22:11:27 ] PROGRAM vipper PARAMETER (n = 5) REAL x(n), y(n) DATA x /0.1 , 0.2 , 0.3 , 0.4 , 0.5 / DATA y /1.228, 1.005, 0.823, 0.674, 0.552/ sumx = 0.0 sumy = 0.0 sumxx = 0.0 sumxy = 0.0 sum = REAL(n) DO 10 i = 1, n sumx = sumx + x(i) sumy = sumy + LOG(y(i)) sumxx = sumxx + x(i) * x(i) sumxy = sumxy + x(i) * LOG(y(i)) 10 CONTINUE det = sum * sumxx - sumx * sumx alog = ( sumxx * sumy - sumx * sumxy ) / det b = ( sum * sumxy - sumx * sumy ) / det WRITE(6, *) ' a = ', EXP(alog), ' b = ', b C check WRITE(6, *) ' i x y', 1 ' a * EXP(bx) (dy)^2' DO 20 i = 1, n yy = EXP(alog + b * x(i)) WRITE(6, *) i, x(i), y(i), yy, (y(i) - yy)**2 20 CONTINUE STOP END >>109　インデント 110 名前：デフォルトの名無しさん mailto:sage [2006/02/05(日) 22:12:16 ] a = 1.499271 b = -1.998701 i x y a * EXP(bx) (dy)^2 1 0.1000000 1.228000 1.227659 1.1664589E-07 2 0.2000000 1.005000 1.005252 6.3688631E-08 3 0.3000000 0.8230000 0.8231379 1.9023346E-08 4 0.4000000 0.6740000 0.6740159 2.5326941E-10 5 0.5000000 0.5520000 0.5519093 8.2298044E-09 Press any key to continue >>110　インデント 111 名前：デフォルトの名無しさん mailto:sage [2006/02/05(日) 22:28:29 ] >>94 死ぬほど易しく作ったから理解しろｗ >>101の本の公式通りに作り、変数名も大体あわせておいた。 ここで誤差のσは各要素で共通一定と仮定した。この場合キャンセルしあうから式から無視してよし。 本を見れば分かるが、二乗誤差が極値をとる条件を、a,bに関してそれぞれ微分を取ることで求めている。 二乗誤差は極値を１個だけもちかつそれは最小だからこれでおｋとなる。 （誤差は外せばいくらでも大きく出来るが、合わせるほうはに限界がある。最高に合いまくって０。） プログラム解説： １０番LOOPで総和を取る。 その後、２＊２の逆行列をかけるのと同じ理屈で、２元連立方程式を解いている。 ２０番LOOPでどのくらい近似できたか見ている。 求めていないが、一番右側の列の二乗差の総和が最小になっているはずである。 a,bの値を手動で変えて、求めた二乗差の総和が最小であることを確かめてみそ。 その結果を考察に書いておけば＋１０点の花丸ハンバーグだ！ 112 名前：94 mailto:sage [2006/02/06(月) 00:15:35 ] みなさんほんと詳しく教えていただきありがとうございます。 ほんと助かりました。 113 名前：デフォルトの名無しさん [2006/02/06(月) 12:36:36 ] ところで、Fortranにオブジェクトデータ構造はいつ入るんですか． 114 名前：sage [2006/02/06(月) 13:32:23 ] >>113 FORTRAN2003。コンパイラーの実装待ち。 115 名前：デフォルトの名無しさん mailto:sage [2006/02/06(月) 17:44:37 ] >>113 オブジェクトデータ構造ってなんですか？ オブイェークト 116 名前：デフォルトの名無しさん [2006/02/06(月) 19:53:17 ] ERK=ER(K)*DCOS(GANMA2(K))-EI(K)*DSIN(GANMA2(K)) EIK=EI(K)*DCOS(GANMA2(K))+ER(K)*DSIN(GANMA2(K)) 実数 ERK=ERK*DCOS(GANMA3(K))-EIK*DSIN(GANMA3(K)) 虚数 EIK=EIK*DCOS(GANMA3(K))+ERK*DSIN(GANMA3(K)) ERK=ERK*DCOS(GANMA4(K))-EIK*DSIN(GANMA4(K)) EIK=EIK*DCOS(GANMA4(K))+ERK*DSIN(GANMA4(K)) この式をEXPを用いた複素数型に変換できなくこまってます。 どなたか優しい方おしえていただけませんか（＞＜） 117 名前：デフォルトの名無しさん mailto:sage [2006/02/06(月) 21:53:43 ] >>116 やぁ！やさしいお兄さんだよ！援助してあげるから、今日からパパと呼ぶんだよ！ 定義：　 ze = er + i ei zg = COS(Γ) + i SIN(Γ) = EXP(iΓ) 積： ze * zg = ( er + i ei ) * ( COS(Γ) + i SIN(Γ) ) 　　　　　= ( er * COS(Γ) - ei * SIN(Γ) ) + i ( ei * COS(Γ) + er * SIN(Γ) ) これより COMPLEX*8 :: ze, zek ze = (er, ei) zek = ze * EXP( CMPLX(0.0d0, gamma2) ) これでＯＫさ！ ただ細かいことを言うと、倍精度複素数はＦ７７標準規格に存在しないので、 完全に規格に乗っ取ろうとすると複素変数は使えないのさ！ でも、大概のコンパイラは倍精度複素数があるのでＯＫさ。 ＥＸＰとＣＭＰＬＸは総称名なので、精度は引数に合わせてくれるはず。 でもマニュアルを調べて総称名が、ちゃんと倍精度複素数になってるかどうか 確認したほうがいいかもね。 118 名前：117 mailto:sage [2006/02/06(月) 22:00:01 ] おっと ze = CMPLX(er, ei) だったね。 さぁ、>>116　おっぱいうｐうｐ 　　　 _ 　∩ 　　(　ﾟ∀ﾟ)彡　おっぱい!おっぱい! 　　(　 ⊂彡 　 　|　　　|　 　 　し ⌒Ｊ 119 名前：デフォルトの名無しさん [2006/02/06(月) 22:39:00 ] なんだかしりませんが Fortran2003って普通にかな〜りすごそうなんですけど ほんとに実装されるんですか？ 120 名前：デフォルトの名無しさん mailto:sage [2006/02/06(月) 23:17:41 ] >>118 つ ttp://kahouweb.hp.infoseek.co.jp/foods/food_230.html 121 名前：デフォルトの名無しさん mailto:sage [2006/02/06(月) 23:29:50 ] >>119 FORTRANの場合スーパーコンピュータやCPUのキラーアプリとして 単独での採算を度外視されて開発されているので、多分実装されるでしょう。 INTELが旧DECのFORTRANコンパイラーチームをCompaqから買い取ったのは そういう方向からだし。 122 名前：デフォルトの名無しさん mailto:sage [2006/02/07(火) 11:50:47 ] 次の宿題もってこい！ 123 名前：デフォルトの名無しさん [2006/02/07(火) 15:28:52 ] ｆ９０を使い初めてプログラムの複素数化をしているのですがうまくいきません、どなたか助言をおねがいします。 PROGRAM TTT IMPLICIT NONE INTEGER :: I,N,M REAL(8) :: EI,ER(0:2047),EI(0:2047),A,B REAL(8) :: FACTR,H,PI,RN1,RN2,EP0,C,WL,TINTY,TH,S(0:2047),TT(0:2047),GEFF,ERI,EII ! PI=4.0*DATAN(DBLE(1.0)) C=3.0E+8 WL=1.55E-6 DO I=0,N-1 FACTR=-H*PI*RN2*RN1*EP0*C/WL TINTY=(ABS(S(I)*S(I))+ABS(TT(I)*TT(I)))/DBLE(2.0) TH=FACTR*TINTY ERI=EXP(-GEFF*H)*(1*DCOS(TH)-1*DSIN(TH)) EII=EXP(-GEFF*H)*(1*DCOS(TH)+1*DSIN(TH)) ER(I)=ERI EI(I)=EII OPEN(20,FILE='TH.dat') WRITE(20,*)ER(I),EI(I) END DO END PROGRAM TTT 上のプログラムからSINやCOSをけしてEXP等で表現 し複素数でのプログラムにするにはどうしたらいいでしょうか？ 124 名前：デフォルトの名無しさん mailto:sage [2006/02/07(火) 16:48:19 ] >>123 まずその前に・・・・・ C=3.0E+8 WL=1.55E-6 これだと単制度で定数が入るので、ここで数値が腐ってしまう。E→Dにせよ。 吸収係数かなんかの計算か？ 発想が逆さまだ。もともとの公式を素直にプログラムするほうがよっぽど楽だぞｗ ほいで (1 + i) * EXP(iθ) = (1 + i) * (cosθ+ isinθ) 　　　　　　　　　　= (cosθ+ isinθ) + (icosθ - sinθ) 　　　　　　　　　　= (cosθ - sinθ) + i(cosθ +sinθ) ゆえに ZE = EXP(-geff * h) * (1.0d0, 1.0d0) * EXP(CMPLX(0.0d0, th)) こうじゃねーか? 間違っていてもシラネｗ 125 名前：デフォルトの名無しさん [2006/02/07(火) 17:30:02 ] >>124 サンクスです、再度考えてみます 126 名前：124 mailto:sage [2006/02/07(火) 17:41:11 ] 補足 (1+i) = SQRT(i) * SQRT(2) = EXP( iπ / 2 ) (i + i) * EXP(iθ) = EXP( i(θ + π/2) ) 127 名前：124 mailto:sage [2006/02/07(火) 17:44:09 ] ごめｗ　今借り物でFEPが違くて勝手が違うｗ　間違い膜リング (1+i) = SQRT(i) * SQRT(2) = EXP( iπ / 2 )　* SQRT(2) (1 + i) * EXP(iθ) = SQRT(2) * EXP( i(θ + π/2) ) つーかやっぱ、一度元の公式に戻れｗ 128 名前：124 mailto:sage [2006/02/07(火) 17:45:34 ] 糞ごめｗ　π/2 →　π/4 うっはーｗｗｗｗｗ　信用度０ｗｗｗｗ 129 名前：124 mailto:sage [2006/02/07(火) 17:52:55 ] ZE = EXP(-geff * h) * (1.0d0, 1.0d0) * EXP(CMPLX(0.0d0, th)) 　 = EXP(CMPLX(-geff * h, th + pi / 4.0d0) こうかな？ｗ スレ汚し　すまんこﾉｼ！ 130 名前：123 [2006/02/07(火) 17:57:00 ] >>124 できましたありがとうございますm(__)m 131 名前：デフォルトの名無しさん mailto:sage [2006/02/07(火) 18:00:55 ] さぁ、>>123　おっぱいうｐうｐ 　　　 _ 　∩ 　　(　ﾟ∀ﾟ)彡　おっぱい!おっぱい! 　　(　 ⊂彡 　 　|　　　|　 　 　し ⌒Ｊ 132 名前：デフォルトの名無しさん mailto:sage [2006/02/07(火) 18:11:54 ] >>131 つ ttp://members.at.infoseek.co.jp/rimssecret/oppai2.htm 133 名前：123 [2006/02/07(火) 18:26:33 ] 実数部分 ERK=ER(K)*DCOS(GANMA2(K))-EI(K)*DSIN(GANMA2(K)) 虚数部分 EIK=EI(K)*DCOS(GANMA2(K))+ER(K)*DSIN(GANMA2(K)) ERK=ERK*DCOS(GANMA3(K))-EIK*DSIN(GANMA3(K)) EIK=EIK*DCOS(GANMA3(K))+ERK*DSIN(GANMA3(K)) ERK=ERK*DCOS(GANMA4(K))-EIK*DSIN(GANMA4(K)) EIK=EIK*DCOS(GANMA4(K))+ERK*DSIN(GANMA4(K)) を同じ考えで EK=E(K)*EXP(CMPLX(0.0d0,(GANMA2(K)+GANMA3(K)+GANMA4(K)))) としたら同結果にならないのですがご教授おねがいします 本当に連投すみませんm(__)m 134 名前：デフォルトの名無しさん mailto:sage [2006/02/07(火) 23:14:45 ] >>133 EK = E(K) * (1.0d0, 1.0d0)**3 * EXP( CMPLX(0.0d0, GAMMA2 + GAMMA3 + GAMMA4) ) = E(K) * SQRT(8.0d0) * EXP(CMPLX(0.0d0, GAMMA2 + GAMMA3 + GAMMA4 + 3.0d0 * pi / 4.0d0) ) でね？ 135 名前：デフォルトの名無しさん mailto:sage [2006/02/08(水) 01:32:58 ] >>133はもう一度>>124を読み直しなさい！！ｗ 136 名前：デフォルトの名無しさん mailto:sage [2006/02/08(水) 13:12:44 ] 誰かネタを投下しろです！ 137 名前：123 [2006/02/08(水) 16:59:13 ] 前半の距離h/2にわたる線形空間伝搬 [exp( j (r2 + r3 + r4 ) h/2 }F{Φ(z,t)}] 　 元の式がこれなのですが複素数でやろうとするとうまくいかないです・・・ (1.0d0, 1.0d0)の部分に変数をいれればいいのですかね？ 138 名前：デフォルトの名無しさん mailto:sage [2006/02/08(水) 17:07:52 ] >>137 その式のままプログラムすればいい。 phi = phi_func(z, t) zf = F(phi) zres = EXP( CMPLX（0.0d0, 0.5d0 * h * (r2 + r3 + r4) ) ) * zf 3行に分けたのはdebug時に、１行ごとにエラーを知るため。 うまく行くのがはっきりしたら、１行にまとめればいい。 139 名前：デフォルトの名無しさん mailto:sage [2006/02/08(水) 18:12:41 ] >>138 コンパイラによってはｆｕｎｃｔｉｏｎを別のfunctionやサブルーチンの引数に取れないかもな。 140 名前：デフォルトの名無しさん mailto:sage [2006/02/08(水) 22:38:28 ] 1 名前：ぶつわよ！φ ★[] 投稿日：2006/02/08(水) 18:26:54 ID:???0 毎日新聞「開かれた新聞」委員会の１２、１月度見解を報告します。 今回は、昨年１２月２５日に山形県庄内町で起きたＪＲ羽越線の特急脱線転覆事故を扱った 毎日新聞の社説「安全管理で浮ついてないか」（１２月２７日朝刊）が 「風の息づかいを感じていれば、事前に気配を感じていたはずだ」と指摘したことなどに対し、 読者から多数の批判意見が寄せられた問題を中心に取り上げます。 ※とっても長いので全文はソースで www.mainichi-msn.co.jp/shakai/wadai/news/20060207ddm012070089000c.html 141 名前：デフォルトの名無しさん mailto:sage [2006/02/09(木) 01:12:17 ] プログラミングの宿題です。助けて下さい。 2種の生物x(t)とy(t)が同一場所に生息している。 x(t)の出生率はa>0であり、x(t)はy(t)に食べられる、その死亡率はby(t)<0であるとする。 一方x(t)を食べるy(t)の死亡率はc<0であり、食料となるx(t)に比例して出生率がdx(t)>0であるとする。 このとき次の連立微分方程式系のシミュレーションのプログラムを作れ。 dx/dt=(a+by)x dy/dt=(c+dx)y ただし、t=0、x=3、y=1、h=0.2、a=1、b=-1、c=-1、d=1(hは時間tの刻み幅) 　 Fortran77の形式でお願いします。よろしくお願いします。 142 名前：デフォルトの名無しさん [2006/02/09(木) 02:18:09 ] >>141 ほれ！ティンティンうｐねｗ 　微分方程式の積分を最高簡単にしておいた。誤差が溜まって結果怪しい。 　ルンゲ＝クッタとか真面目なのに自分で治しとけｗ PROGRAM oppai PARAMETER (h = 0.2, a = 1.0, b = -1.0, c = -1.0, d = 1.0) t = 0.0 x = 3.0 y = 1.0 DO 10 i = 1, 100 t = i * h dx = (a + b * y) * x * h dy = (c + d * x) * y * h x = MAX(0.0, x + dx) y = MAX(0.0, y + dy) WRITE(6, *) x, y, t, i 10 CONTINUE END >>142　インデント 結果はというと、ｔ＝１０くらいで、生命体ｘが滅亡し、その後餌のなくなった 生命体ｙもすぐに滅亡したｗ 143 名前：142 mailto:sage [2006/02/09(木) 02:27:24 ] >>141 メッシュを細かくして計算したら、ｘとｙは交互に増えたり減ったりして振動したｗ この結果がもっともなので、>>142のままではやはりﾔｳﾞｧｽｗｗｗｗｗｗｗｗｗｗ 144 名前：コンパイル？ [2006/02/09(木) 05:24:58 ] fortran(パソコン？）に初めて触ると言うくらいの初心者です。 www.coastal-env.k.u-tokyo.ac.jp/koibuchi/fortran95/fortran77.htmにそって、".Net"が無い為、fortran77をダウンロードしました。 書いてあるとおり、"F77.exe"の位置を指定して、プログラムをコンパイル(実行）を行った結果、画面には… 「コンパイルできません。 ファイル" …(ファイル名).exe"は、存在しません。」 と表示され、何も反応してくれません。 どうしたらよいでしょうか？教えてください…。 145 名前：コンパイル？ [2006/02/09(木) 05:50:14 ] ちなみに、f77.exeは、"C:\Program Files\Fortran77\ftn77.exe"に保存してあります。 保存先は、"…\My Documents\Fortran"に指定しました。 F95でも、同様な事をすれば"fcpad231.lzh"でうまく行くと書いてありますが、こちらも教えて頂けましたら、幸い名のですが…。 (今、自分は研究者なのですが…遅まきながらの数値計算の導入で、かなり四苦八苦しています。よろしくお願いします。） 146 名前：デフォルトの名無しさん mailto:sage [2006/02/09(木) 07:03:53 ] >>144 それは 「あなたの作ったプログラムに間違いがあるのでコンパイルできませんでした。 （コンパイルできなかったので）実行ファイルが存在せず、実行できません」 という意味です。 「メッセージ」の欄にエラーの起こった行とその理由が表示されるので それを参考にもう一度よくプログラムを見直して訂正して下さい。 >>145 FTN95 コンパイル時のオプション指定の方法が FTN77 と違うため、 そのままでは CPad for Salford FTN77 では使えません。 FTN95付属のIDEを使うか（但し、日本語は使えません）、 適当なエディタを使って下さい。 147 名前：コンパイル？ [2006/02/09(木) 07:52:35 ] >>146 ありがとうございます。そうだったんですね。 メッセージには、 ■C:\Documents and Settings…\デスクトップ> ftn77 /link ex.for [CPad: コンパイル終了] とあります。 どうしたらよいのでしょうか…？ あと、適当なエディタを用いるとありますが、どのようなものを使用すればよいのでしょうか？ 本当に初歩的な質問ですいません…。 148 名前：デフォルトの名無しさん mailto:sage [2006/02/09(木) 08:03:53 ] もう一つ質問させてください。 f77よりf90でパンコイルしたほうが高速ですか？ 149 名前：デフォルトの名無しさん [2006/02/09(木) 11:33:46 ] >>148 初心者レベルで問題になるような著しい差はない。 エディターに関しては、初心者だというなら 何がしたいのか分からないが、新たに言語を習得したいというならF90を 学ぶほうがいい。 150 名前：デフォルトの名無しさん mailto:sage [2006/02/09(木) 11:38:41 ] エディターに関しては、初心者だというなら予約語（？）が色づけされる タイプのものがいいはずなので、統合環境のものを使うのが吉だろう。 MSの開発環境は複雑化しすぎて初心者には荷が重いだろう。 無料のものといえば、旧ソ連の金のない連中が、g77や著作権切れの 本のPDFなどをセットにしたFORTRAN統合環境を出していたがｗｗｗ ま、FTNでいいんでないの？ 151 名前：デフォルトの名無しさん mailto:sage [2006/02/09(木) 20:29:47 ] >>147 >どうしたらよいのでしょうか…？ あなたが使っているPCの環境やプログラムの内容がわかるエスパーが現れるまで しばらくお待ち下さい。 152 名前：デフォルトの名無しさん [2006/02/10(金) 02:28:09 ] >>141　今日は2次のルンゲ=クッタにしてみたｗ PROGRAM oppai2 PARAMETER (h = 0.2, a = 1.0, b = -1.0, c = -1.0, d = 1.0) t = 0.0 x = 3.0 y = 1.0 hh = h / 2.0 DO 10 i = 1, 200 t = i * h dx0 = (a + b * y) * x * hh dy0 = (c + d * x) * y * hh dx = (a + b * (y + dy0) ) * (x + dx0) * hh dy = (c + d * (x + dx0) ) * (y + dy0) * hh x = MAX(0.0, x + dx) y = MAX(0.0, y + dy) WRITE(6, *) x, y, t, i 10 CONTINUE END >>152　インデント とりあえず、振動が続くようになった。捕食のモデルから予想されるものに なっていると思う。 ま、普通は４次のルンゲ=クッタなんだがｗ 153 名前：デフォルトの名無しさん mailto:sage [2006/02/10(金) 02:47:02 ] >>152　4次ルンゲ食った PROGRAM oppai4 PARAMETER (h = 0.2, a = 1.0, b = -1.0, c = -1.0, d = 1.0) t = 0.0 x = 3.0 y = 1.0 hh = h / 2.0 DO 10 i = 1, 2000 t = i * h dx1 = (a + b * (y ) ) * (x ) * hh dy1 = (c + d * (x ) ) * (y ) * hh dx2 = (a + b * (y + dy1) ) * (x + dx1) * hh dy2 = (c + d * (x + dx1) ) * (y + dy1) * hh dx3 = (a + b * (y + dy2) ) * (x + dx2) * hh dy3 = (c + d * (x + dx2) ) * (y + dy2) * hh dx4 = (a + b * (y + dy3) ) * (x + dx3) * hh dy4 = (c + d * (x + dx3) ) * (y + dy3) * hh x = MAX(0.0, x + dx1 / 6.0 + dx2 / 3.0 + dx3 / 3.0 + dx4 / 6.0) y = MAX(0.0, y + dy1 / 6.0 + dy2 / 3.0 + dy3 / 3.0 + dy4 / 6.0) WRITE(9, *) x, y, t, i 10 CONTINUE END >>153　インデント 振動しつつ、ｘ、ｙとも１に収束してゆく。　ｔ＝３５０　くらいでほぼ収束か。 154 名前：デフォルトの名無しさん mailto:sage [2006/02/10(金) 22:57:53 ] ポトペタ可能な開発環境はありませんか？ 155 名前：デフォルトの名無しさん mailto:sage [2006/02/10(金) 23:17:14 ] ポトペタ　ってなんですか？ 156 名前：デフォルトの名無しさん mailto:sage [2006/02/10(金) 23:34:45 ] >>154 こんなんがあるな。使ったこと無いけどｗ ttp://www.hulinks.co.jp/software/ginomenustudio/index.html つーか、GUI は他の言語で作った方が楽だけどな。 157 名前：デフォルトの名無しさん [2006/02/11(土) 14:41:04 ] >>147 ・・・いっぽうロシア人は鉛筆を使った。・・・・ 金なしロシアのｆ７７開発環境だｗ www.imamod.ru/~vab/vfort/ 158 名前：デフォルトの名無しさん [2006/02/11(土) 22:26:04 ] 多変数多項式のGCD計算がなかなか作れないのですが、どなたか見本みたいな感じでソースを見せて頂けませんか？ あと、グレブナ基底を求めるプログラムもあれば頂きたいです。 159 名前：デフォルトの名無しさん mailto:sage [2006/02/11(土) 22:48:29 ] >>158 　【FORTRAN７７】　ﾌﾟﾛｸﾞﾗﾑを起動させて！！ science4.2ch.net/test/read.cgi/sim/1057990786/68 はい、次の方どうぞー 160 名前：デフォルトの名無しさん mailto:sage [2006/02/12(日) 01:44:35 ] >>158　いろいろあった気がするが。　　岩波から本が出てなかったか？　買え！ To search the Netlib repository, enter your search query in the form below and then press the Go button. Look for : any of the wordsall the words Query "groebner" found 1 matches. 1. toms/628 keywords: groebner basis, polynomial ideals, rational integers gams: C3b title: GROEB for: canonical (or Groebner) bases of polynomial ideals by: F. Winkler et al. ref: ACM TOMS 11 (1985) 66-78 Score: 100% www.netlib.org/cgi-bin/search.pl 161 名前：デフォルトの名無しさん [2006/02/12(日) 02:11:34 ] すごい資料が出たらしいぞｗｗｗｗ bbs.enjoykorea.jp/tbbs/read.php?board_id=thistory&nid=1578859 これをみれ！　そして１００回saveしろ！ 162 名前：158 [2006/02/12(日) 02:22:29 ] >>159 あ〜、すみません。それ俺です… 向こうは人が少ないし、偶々ここを見つけて、こっちの方が本職かな？と思いまして。 >>160 岩波からですね。漁ってみます。 自分は脳ミソの出来がよろしくないので、わからないときはまた質問させて頂きます m(_ _)m 163 名前：158 [2006/02/12(日) 02:39:54 ] 連続ですみません。 >>160 この検索すごいですね。grobnerもGCDも一発で見つかりました。 まだ中身を細かく見てないですが、 岩波の本を買って参照しながら解読してみます。 助かりました。 164 名前：デフォルトの名無しさん mailto:sage [2006/02/12(日) 02:55:11 ] >>162 岩波でなかったかもｗ　間違ってたらすまんｗ プログラムはよくネットに落ちているので検索すればもっと出てくると思う。 165 名前：デフォルトの名無しさん [2006/02/12(日) 02:59:24 ] >>163 NetLibはインターネット以前からある数学系のサブルーチン集だ。 古いのがごみのように詰まっているｗ　 166 名前：158 [2006/02/12(日) 04:18:13 ] >>164 それらしい本で 岩波から計算代数と計算幾何(佐々木,今井,他,共著) 共立から計算機数学(野崎著) を本屋の検索で発見したので、明日図書館で調べてみようかなと思います groebnerのプログラムが長すぎて脳がスーパーノヴァしてるので、今日は寝ることにします(_ _;) 数値解析の本ばかりで、GCD計算は代数系の入門書からだったので、本当に助かりました >>165 古くても重宝しますね、これは。 数にビビリました。 ルンゲも、ガウスも、DE公式も、SVDも、なんでもありますね！ 無償で頂けるなんて親切な人達がいるもんですね(^-^) 167 名前：デフォルトの名無しさん mailto:sage [2006/02/12(日) 17:15:33 ] >>166 昔の大型計算センターはどこでもライブラリが大量に公開されているのが普通だった。 ＦＯＲＴＲＡＮの世界は歴史が古いので、ソフト囲い込み運動と その反動のオープンソース運動とは無関係にソース公開共有の 伝統がまだ残っている。田舎で鍵をかけないで出かける世界？ｗ 168 名前：デフォルトの名無しさん [2006/02/14(火) 03:02:05 ] さぁ次の宿題もってこい！ 169 名前：デフォルトの名無しさん [2006/02/14(火) 20:32:01 ] 配列を返す関数はどのように定義するのですか？ 170 名前：デフォルトの名無しさん mailto:sage [2006/02/14(火) 22:20:04 ] >>169 program hoge 　implicit none 　real, dimension(3) :: x, y, z 　real, dimension(5) :: u, v, w 　x = (/1, 2, 3/); y = (/2, 3, 4/) 　u = (/4, 5, 6, 7, 8/); v = (/10, 20, 30, 40, 50/) 　z = adda(x, y) 　w = adda(u, v) 　print *, z 　print *, w 　stop contains 　function adda(a, b) 　　real, intent(in), dimension(:) :: a, b 　　real, dimension(size(a)) :: adda 　　integer :: na, nb 　　na = size(a); nb = size(b) 　　if (na > nb) then 　　　adda(1:nb) = a(1:nb) + b(1:nb); adda(nb+1:na) = 0.0 　　else if (na < nb) then 　　　adda = a(1:na) + b(1:na) 　　else 　　　adda = a + b 　　end if 　end function adda end program hoge 171 名前：169 mailto:sage [2006/02/14(火) 23:25:44 ] >>170 レスありがとうございます！わかりました！ 172 名前：デフォルトの名無しさん mailto:sage [2006/02/15(水) 19:35:50 ] pathとファイル名指定して、そこでファイルオープンするような処理を 書きたいのですけど、char型に突っ込んだ後の、うしろの余ったスペースが どうにもならないっす。 みなさんはうしろの余ったスペースどうやって消してますか？ 173 名前：172 [2006/02/15(水) 20:14:13 ] age忘れたので、ageさせてください。すいません。 174 名前：デフォルトの名無しさん mailto:sage [2006/02/15(水) 20:27:50 ] >>172 ! TRIM 関数, LEN_TRIM 関数の使い方 program huge implicit none character(len=30) :: path, fname character(len=70) :: r1, r2 print *, "path:"; read *, path print *, "file name:"; read *, fname r1 = path // fname ! そのまま連結 r2 = trim(path) // fname ! 末尾の空白を除いて連結 print *, r1 print *, r2 print *, len(r2) ! 文字列の長さを返す print *, len_trim(r2) ! 末尾の空白を除いた文字列の長さを返す stop end program huge 175 名前：172 mailto:sage [2006/02/15(水) 23:55:23 ] >>174 わざわざプログラムまで書いていただいて、ありがとうございます！ おかげでわかりました！ 176 名前：デフォルトの名無しさん [2006/02/17(金) 09:54:08 ] 1 FORMAT(1H ,I5) 2 FORMAT(1H ,E15.7) 3 FORMAT(1H ,I5,1H ,I5,1H ,7E15.6) 4 FORMAT(1x,I5,1X,I5,1X,E15.6,E15.6) 5 FORMAT(B20) 6 FORMAT(I5,1H ,I5) 7 FORMAT(1H ,I5,1H ,3E15.6) 11 FORMAT(7X,I5,1X,E15.6,E15.6) fortran90からはじめたのですが上のFORMATｎ意味がよくわかりません どなたかご教授おねがいします 177 名前：デフォルトの名無しさん [2006/02/17(金) 10:43:09 ] HはホレリスのH。廃止勧告が出ていて、いまではHな人しか使わないｗ 基本的には文字列を囲む引用符”　”の祖先で、以降が文字列ということを示している。 文字数はHの前の数字で表している。 FORTRAN66までは文字型が無かったという事情がある。 １文字目が１Hになっていることが多いのは、かつてのFORTRANでは出力の１文字目は プリンター制御文字になっているという事情がある。改行、改ページ、改行せずができる。 （ビデオコンソールの普及は７０年代以降なのでコンソール出力を前提とした発想は捨てよ）。 それゆえ、この１文字目を空白で飛ばさないと、頭の１文字が消えてしまうという現象がおきる。 いまでもコンパイラーごとに１文字目の扱いは異なる。VisualFOTRANなどではコンパイラ オプションで１文字目の扱いを変えることができる。 最近のコンパイラーは画面出力では１文字目をプリンター制御文字と解釈しない傾向がある。 そういうコンパイラーで実行すると、常に行頭に１文字空白が付いたり、時々０，１が現れることになる。 178 名前：デフォルトの名無しさん [2006/02/17(金) 10:59:47 ] プリンター制御文字　４通り ０　２行改行 １　改ページ ＋　改行せず その他　１行改行（普通の出力） その他 Xは入力では頭の数だけカラムをスキップ、出力では頭の数だけ空白出力。 Bは２進法で出力。B２０なら２０桁の幅で2進法出力。 I,ESは基本なのでマニュアル嫁。 １H　が冒頭に連発するのは初心者には理解できない躓きポイントではあるが、 これは上代の古語が出てきたと思って、ももしきや〜と古きをしのんで喜ぶべき ポイントでもある。歴史と伝統のFORTRANの醍醐味でもある。 あなとうと、あなうるわし。 179 名前：デフォルトの名無しさん mailto:sage [2006/02/17(金) 11:57:04 ] FORMATは昔からFORTRANの難所の一つとされてきたが、 最近はFORMATのほうも進化してきたので学び甲斐もあるべ。 配列出力で、：や（）による反復指定を使って思い通りに出力できると、 他の言語では得られない快感が得られる。 180 名前：176 [2006/02/17(金) 18:50:07 ] ありがとうございました＾＾ がんばって理解していこうと思います 181 名前：デフォルトの名無しさん [2006/02/19(日) 10:26:10 ] COMMON /M/!,N,LZ COMMON /P/P,GAMMA,GPH とゆうCOMMON文があるのですがCOMMON文自体を理解していなく 何の処理をしているのかいまいちわかりません ｆ９０ではMODULE文にできるらしいのですがどうのように 書き直せばいいのかどなたかご教授お願いします。 182 名前：デフォルトの名無しさん mailto:sage [2006/02/19(日) 12:47:24 ] >>181 COMMONブロック P の先頭から順に実数型変数 P, GAMMA, GPH が割り当てられる。 （同じプログラム単位中でP, GAMMA, GPHの定義がない場合） COMMONブロック M は変数名の並びがコメントアウトされているので無効？ （エラーになる？） モジュールにするとこんな感じ module hoge implicit none private !integer, public, save :: N, LZ real, public, save :: P, GAMMA, GPH end module hoge 但し、他のプログラム単位中でCOMMONブロックM, Pが他の用途に使われていたり P, GAMMA, GPHの並びが違っていたりすると正常に動作しないことがある。 183 名前：デフォルトの名無しさん [2006/02/19(日) 16:57:35 ] COMMON文はグローバル変数としての機能もあるが、その概念を外延すると 理解しにくいこともある。むしろ同じメモリー領域の手動割付と考えると分かりやすい。 昔メモリーが少なくて貴重だった時代には、手動で半動的にメモリー割付を してメモリーを有効活用していた。こういう場合同じCOMMONブロックに全然 対応の付かない変数が書かれていたりする。 COMMONブロックでは、書いた変数の名前の順に先頭から順番にメモリーが 割り振られてゆく。 COMMON /P/P,GAMMA,GPH この場合Pという名前のブロックに、（P,GAMMA,GPHを単精度４バイトとして） P,GAMMA,GPHに対応する４バイトの領域が３つとられる。 このコモンブロックを別のサブルーチンでは、COMMON/P/A,B,C　と受けることも 出来る。この場合P,GAMMA,GPHの内容はA,B,Cと一致する。（変数の番地が 一緒になる） またこのコモンブロックを別の形で受けてもいい。 COMMON/P/A(3) で受けた場合、P,GAMMA,GPHの内容はA(1),A(2),A(3)に対応する。 重要なのは順番で変数名ではない。これはMODULEの場合と異なる。 COMMON /P/GAMMA,GPH,P で受けたとすると、元のP,GAMMA,GPHはそれぞれGAMMA,GPH,Pに対応する。 このように違う名前でCOMMMONを受けていなければ、単純にMODULEに直せる。 そうでない場合は慎重に作業する必要がある。 184 名前：デフォルトの名無しさん [2006/02/19(日) 23:37:30 ] この式のプログラムできませんか？さっぱりわからんのです。。。 F(P)=(9/(128*PAI**2))*(1/P)|EXP(-P*(1/TAN(X/2)))*(1+COSX**4)dX [|]はインテグラルで、範囲は０からπ/２です。シンプソンのプログラムはあるんですが、もうさっぱりです！ どなたかお願いします。。。 185 名前：デフォルトの名無しさん [2006/02/20(月) 02:13:21 ] >>184　 tan(X/2)=1/Yとおくと積分範囲は∞〜１ dX=-2/(1+Y^2)dY cos(X)=(Y^2-1)/(Y^2+1) 1/(1+cos(X)^4) = (Y^2+1)^4/(Y^8+6Y^4+1) 代入して (1/p)∫exp(-p/tan(X/2)) / (1+cos(X)^4) dX = -(2/p)∫(Y^2+1)^3 exp(-pY) / (Y^8+6Y^4+1) dY 積分範囲を１〜∞に反転させれば頭のマイナスが消える。 とりあえず原点付近の１／ｔａｎの発散を避けるようにしてみた積分範囲が無限大にまでになるが、 被積分項は急速に小さくなるので適当なところで打ち切れるだろう。 適当なので間違ってたらスマソ。 186 名前：デフォルトの名無しさん [2006/02/20(月) 12:18:04 ] fortranは、コマンドライン引数使って実行できないのですか？ 187 名前：デフォルトの名無しさん mailto:sage [2006/02/20(月) 12:23:48 ] >>186 iargc() getarg 188 名前：デフォルトの名無しさん mailto:sage [2006/02/20(月) 12:35:30 ] >>186 規格に入ったのはFortran2003から、Fortran95までの規格には存在しない。 しかしたいていはベンダー拡張ルーチンでできることが多い。 マニュアルを調べてみるといい。POSIX互換の形式をとることが多いので、>>187が 書いたようなコマンド名になることが多い。 189 名前：デフォルトの名無しさん mailto:sage [2006/02/20(月) 15:16:01 ] >>187, 188 なるほど、勉強になります。ありがとうございます。 190 名前：デフォルトの名無しさん [2006/02/20(月) 18:40:38 ] >>182 183 ありがとございます m(__)m 191 名前：デフォルトの名無しさん [2006/02/20(月) 23:15:53 ] VAXVMS上で動く，Fortran IVで書かれたプログラムを Fortran77以降のものに自動変換してくれるプログラムとか 無いですか？ 192 名前：デフォルトの名無しさん mailto:sage [2006/02/21(火) 00:07:10 ] >>191 見たこと無い。 193 名前：デフォルトの名無しさん mailto:sage [2006/02/22(水) 13:44:48 ] 宿題を頼んだ人は事後報告をするように。 194 名前：デフォルトの名無しさん mailto:sage [2006/02/22(水) 19:59:20 ] >>193 宿題に困るレベルの人が総括（事後報告）できると思う？ 「ありがとうございます」で納得してあげてちょ。 195 名前：デフォルトの名無しさん mailto:sage [2006/02/23(木) 00:17:12 ] 解答もらったらそれっきり音沙汰ない奴もおるからのぉ。 196 名前：デフォルトの名無しさん [2006/02/24(金) 01:07:22 ] さて、遠慮せず次の宿題を持ってきたまへ。 197 名前：デフォルトの名無しさん mailto:sage [2006/02/24(金) 01:31:52 ] 宿題たまってるけど眠いのでもう寝ます。 お休みなさい。 198 名前：デフォルトの名無しさん [2006/02/24(金) 07:10:10 ] Fortranで書かれた数値計算のプログラムをC#に移植しようとしています。 至る所で引数に配列の要素数をもつ関数が定義されているのですが、 C#のGetLengthやLengthをつかってその都度判定するような仕様にしては問題あるのでしょうか。 つまり、GetLengthなどの関数を実行する時間を節約するために要素数を引数で与えているのでしょうか。 それともFortranでは配列数を取得するような関数が無いためなのでしょうか。 199 名前：デフォルトの名無しさん [2006/02/24(金) 12:39:09 ] >>198 FORTRNでは、FORTRAN９０になるまで配列の大きさを知る命令は存在しない。 FORTRAN９０で関数SIZEが導入された。 FORTRAN７７以前では、配列の大きさは副プログラム（サブルーチン、関数）に、 引数として明示的に渡すか、COMMON等で渡すか、あるいは”知っている”もの としてサブルーチン側で勝手に使うかした。 FORTRANでは変数はアドレス渡しであり、FORTRAN７７以前の書式を使うと 配列であっても渡しているのは先頭要素のアドレスだけで、配列の大きさもの情報も 次元も渡していない。（FORTRAN９０で導入された新しい書式を使うと配列の サイズや次元の情報も渡される。） ここで気をつけなければならないのは、FORTRAN６６時代のプログラムでは、 サブルーチン側の引数の仮配列の大きさを実際の長さとは無関係に１０などにして 内部では好き勝手な長さとして使用する書き方が標準的に採用されていたことだ。 （FORTRAN77では、サイズを＊で定義するやり方に相当するものかな？） こういう書法の伝統を引いている場合、>>198の移植方法では大混乱をきたすだろう。 200 名前：199 [2006/02/24(金) 12:53:08 ] 参考　（ｔabが消えるので７カラム目までの空白は脳内で補完してください） ＊FORTRAN66時代にありがちなプログラム REAL A(20, 20), B(200) CALL ZERO(200, B) CALL ZERO(400, A) CALL ZERO(20, A(1, 2)) CALL ZERO(20, A(1, 3)) CALL ZERO(20, A(1, 4)) END SUBROUTINE ZERO(N, X) REAL X(10) DO 10 I = 1, N X(I) = 0.0 10 CONTINUE END これは昔よく使われたやり方。サブルーチンに渡されるのは配列の先頭要素の アドレスだけなので、次元が整合していなくても問題はない。また配列の途中 だけを渡して部分的に利用しても問題はない。（FORTRANの配列はCとは 逆順で列方向に対して連続的） CALL ZERO(200, B)　これは素直な呼び方。 CALL ZERO(400, A)　これは２次元配列を１次元的に還元して渡している。 CALL ZERO(20, A(1, 2))　これは２次元配列の２列目だけを渡している。 CALL ZERO(20, A(1, 3)) この呼び方は列ベクトルの使用法として良く使われる。 CALL ZERO(20, A(1, 4))　同様に３，４列目を渡している。 ＊FORTRAN77時代 SUBROUTINE ZERO(N, X) REAL X(N) または REAL X(*) 201 名前：199 mailto:sage [2006/02/24(金) 13:05:48 ] １９９，２００で書いたこととは無関係に一般論的な結論を言わせてもらうが、 fortranでは宣言した配列の大きさとは無関係な大きさで配列を使うことが ままあるので、引数として渡している大きさを配列の宣言長に直すのは、 とても危険である。 FORTRANのライブラリなどでやたらと配列の寸法を引数として求められるのも、 配列は先頭要素のアドレスしか渡しておらず、そのほかの情報を一切渡していない からである。 またFORTRANでは、この性質を大いにプラスに利用しているので （というかこれが無いとまともな実用的なプログラムが書けない）これを非難するのは お門違いである。 FORTRAN９０では、配列の次元や寸法も渡しているがその性で、>>200で書いた ような列ベクトルを渡す時に、余分な処理が必要になって速度が低下する。 こういう時にFORTRAN77式の先頭アドレス渡しの簡潔さが心地よく感じられる。

---

---

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

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