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

くだすれFORTRAN（超� ..

620:617
06/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
06/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:デフォルトの名無しさん
06/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:さおり
06/07/12 16:17:47
すいません！無理を承知でカキコします。
私は今理系で某国立大学に通ってるんですけど、水商売のバイトにハマってしまって全然授業についていってない状態なんです。
それでプログラミングの授業があるんですけど課題が出て「放物の運動（空気抵抗を考慮する）」のプログラムをfortranで作れって課題が出たんです。。
初速度と投げる角度を入力して一定時間ごと（１秒ごととか）の物体の位置をｘとｙで出すらしいんです。
私、授業にほとんど出てなくて本当にまったくわかってない状態なんです。
少し真面目な大学なんで私みたいな子は浮いてしまって、学校には助けてくれる友達もいないんです。
誰かプログラム作ってください。お願いします。この授業を落としてしまうと本当に留年なんです。おねがいします！！

624:デフォルトの名無しさん
06/07/12 16:43:42
>>623
まずはおっぱいうｐ汁

625:デフォルトの名無しさん
06/07/12 21:23:23
>>623
担当教官をお店で接待すれば優もらえるおｗ

626:デフォルトの名無しさん
06/07/13 00:14:44
金があるなら留年しろ。

627:デフォルトの名無しさん
06/07/13 01:04:04
俺は少年なら助けてやるのだが・・・・女は興味ないｗｗｗ

628:デフォルトの名無しさん
06/07/13 07:18:28
>>623
課題だけ書く方が答えが得られるすれのようです。

629:さおり
06/07/13 10:47:04
本当にお願いします！
いろいろ事情があって留年できないんです。

630:デフォルトの名無しさん
06/07/13 12:22:06
>>629
課題は省略しないで出題されたまま提示するべき
ヒントになるかどうか知らん
URLﾘﾝｸ(www12.plala.or.jp)

631:デフォルトの名無しさん
06/07/13 12:26:25
>>623
空気定数って物体のサイズによって二通りの計算方法があるけど、どっちを使うの？

632:631
06/07/13 12:29:35
>>629
>>630 のヒントの方の計算方法だと楽でいいよね。さあどっちだろうか？

633:631
06/07/13 16:25:24
>>631
×空気定数→○空気抵抗

634:さおり
06/07/13 21:43:35
えーっと空気抵抗が速度の２乗に比例する場合らしいです。本当にありがとうございます。

635:デフォルトの名無しさん
06/07/13 22:08:31
解決したみたいだから、次の宿題。

636:デフォルトの名無しさん
06/07/13 22:56:22
これから宿題をストーリーつきで頼む人は、

『ボク宿題が解けなくて困ってるんだ。
お兄ちゃん様に一緒に解いてほしいなぁ～。云々』

という、もっとやる気を出させるストーリーを考えること！

637:デフォルトの名無しさん
06/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:デフォルトの名無しさん
06/07/13 23:19:43
>>636
兄貴ぃ～！
俺、宿題分かんないッス！！
徹夜で手取り足取り教えて欲しいッス！！

639:デフォルトの名無しさん
06/07/14 00:15:18
>>638
サブ！　手取り腰取り教えてやるぜ！！

640:デフォルトの名無しさん
06/07/14 01:05:54
>>637
最後まで解けよｗ

641:デフォルトの名無しさん
06/07/14 02:10:46
>>637
言いだしっぺのh(ry

642:641
06/07/14 02:11:58
間違った
637でなくて>>640

643:637
06/07/14 10:13:37
>>640
そうじゃなくて、学校で習った数値解析法で解かないと
合格になる気がしないんですけど。XX法とか言ってもらえれば
続きも出来るのですが……

644:デフォルトの名無しさん
06/07/14 11:08:41
>>634の書き込みを見ると、「らしい」なんて書いているから、課題の正確な内容を
把握してなさそう。仮に折角>>637が解いたところで「やっぱり違った」とかで別解を
要求してきそうだ。

645:デフォルトの名無しさん
06/07/14 13:55:36
C++スレ見たいに宿題依頼のテンプレが必要なんでね？

見たところ、宿題はいまだF77中心のようだし。ｆ９０かｆ７７かも選択してもらわないと。

646:さおり
06/07/14 16:17:43
そんなこと絶対しません。よろしくおねがいします。

647:デフォルトの名無しさん
06/07/14 18:27:20
つーことで、さおりタンの仕事は
1. 微分方程式の確認
2. 授業でならった微分方程式の数値解析法の(せめて)名前の確認
だね。

648:デフォルトの名無しさん
06/07/14 19:44:38
>>646
「そんなことしません」じゃなくて、レポートとして期待されているプログラムを
作成するためには情報が足りねーの。
課題を入手してきなさい。あと、おっぱいうｐな。

649:デフォルトの名無しさん
06/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
06/07/14 21:17:58
>>649
thx.
いやー。数値解析の教科書とか見て気が重かったのだが、ありがとう。
さおりタン満足してくれるといいね。

651:649
06/07/14 21:40:27
>>650
>>649はあまりに手抜きｗ　誰かもっと良いの作ってやれｗ

２元２階の微分方程式なので、本のルンゲ・クッタルーチンにそのまま放り込めないｗ

ちょっと脳みそ使う必要がある。数独１回分くらいｗ　

今日暑いから考えらんねー　冷房ほしいよー

652:デフォルトの名無しさん
06/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:デフォルトの名無しさん
06/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:６５３
06/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:デフォルトの名無しさん
06/07/15 03:04:03
　　　　　 |　|/ノ二__‐─ｧ　　ヽﾆニ二二二ヾ } ,'⌒ヽ
　　　　 /⌒!|　 =彳ｏ｡ﾄ￣ヽ　　　　　'´ !o_ｼ｀ヾ　| i/ ヽ !　　　　ウホッ！　いい男・・・
　　　　 ! ハ!|　　ー─　'　　i　　!　　 `' 　 '' " 　 ||ヽ　l |
　　　　| | /ヽ!　　　　　　　　|　　　　　　　　　　　 |ヽ i　!
　

656:652
06/07/15 09:43:09
>>653
ウホッ！

657:デフォルトの名無しさん
06/07/15 09:44:53
>>623,>>629,>>634,>>646

>>637をもとに>>649が解いたものは、m、k、初期位置、時間の刻みに
適当な値を代入して解いている。本来はあんたが提示すべき情報。

>>630,>>647,>>648を良く読め。

658:デフォルトの名無しさん
06/07/15 10:57:31
教えて！goo にも同じような質問が・・・

659:637
06/07/15 12:20:32
>>658
じゃあ我々の役目も終ったのかな？
みんなお疲れ様。

660:デフォルトの名無しさん
06/07/15 17:35:26
女が来るとふいんきが乱れるから、ホモOnlyスレにしないか？

661:デフォルトの名無しさん
06/07/15 18:59:43
プログラムの実行中に一瞬ポーズする関数ってありますか？

662:デフォルトの名無しさん
06/07/15 19:17:53
>>661
PAUSE

663:デフォルトの名無しさん
06/07/15 19:30:41
>>662
ただPAUSEはFORTRAN2003あたりで廃止勧告が出されている。
TSS利用で便利だから漏れは好きなんだけどなｗ

664:デフォルトの名無しさん
06/07/15 20:28:29
write(*,*) abs(-16.8-(6.8))>10.0
を出力すると「T」になってしまいます。
助けてください。

665:664
06/07/15 20:32:23
間違えました、
real*8 a,b
a=-16.80000000000000
b=- 6.80000000000000
write(*,*) abs(a-b)>10
とやったら「T」と出てしまいます。
何が原因でしょうか？

666:664
06/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:デフォルトの名無しさん
06/07/15 20:59:25
FORTRAN77は大文字でプログラム書くのがデフォルトなのか

668:デフォルトの名無しさん
06/07/15 21:04:15
>>664
それは良くあることｗ　気にするな。　浮動小数点之罪。

669:デフォルトの名無しさん
06/07/15 21:06:04
>>667
本来のFORTRAN77の規格では大文字のみ。

そもそも、昔は小文字が表示できない端末とか良くあったｗ

670:デフォルトの名無しさん
06/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
06/07/15 21:45:51
>>670
どうもありがとうございます。
一応計算中の変数を画面に表示させてチェックしたのですが、
aは-16.8000000000000
bは -6.8000000000000
と表示されていました。
この上でabs(a-b)>10をやるとTになりました。

672:デフォルトの名無しさん
06/07/15 22:03:21
>>671
(T_T)
30桁くらい出力してごらん。
write ( *, '(2f35.30)' ) a, b

673:670
06/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:デフォルトの名無しさん
06/07/16 14:38:42
　積分区間の下限Aと上限B、分割数N（偶数）を読み込んで、
Simpson公式を用いて、定積分I=(A,B)∫f(x)dxを計算するプログラムを作れ。
　ただし、数値積分の部分のプログラムは、任意の関数に対して使用できる
サブルーチン・サブプログラムを用いよ。その場合関数f(x)の計算は関数サブプログラムで行う。

　さらに上記の分割数Nが奇数に対しても使えるように拡張せよ。
Simpson公式を用いて、Nが奇数のときは1区画だけ台形公式を用いるようにすればよい。

手も足も出ません…、お願いできますか？

675:デフォルトの名無しさん
06/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:デフォルトの名無しさん
06/07/16 21:55:10
ファイルを読み込んだり作ったりするとき、
実行ファイルと同じディレクトリの中に
別のディレクトリを作って保存したりするにはどうしたら良いでしょうか。

677:デフォルトの名無しさん
06/07/16 22:07:35
>>676
既出の質問だが・・・・お答えしよう。

FORTRANの標準命令ではできない。

コンパイラメーカー固有の拡張命令でシステムの関数を呼び出すものがあるだろうから
それらを通じてやるしかない。マニュアルの後ろのほうの章を探して味噌。

678:676
06/07/16 23:32:02
>>677
ｻﾝｸｽｺ

679:デフォルトの名無しさん
06/07/17 12:21:13
番組の途中ですが・・・・・・・・・・・・

盧武鉉（ノ・ムヒョン）大統領が今月11日に行われた与党ヨルリン・ウリ党の指導部
および国会の統一外交通商委員会に所属する議員らとの晩さん懇談会で行った
発言が波紋を呼んでいる。

一部新聞は懇談会出席者の証言を引用し、盧大統領は「ブッシュ米大統領が
北朝鮮問題を善と悪の対立概念で見ているため、説得が難しくなっている。
米国は友邦なので厳しく責めることは出来ないが、
日本とは対決しなければならない」と語ったという。
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

ソース:朝鮮日報
URLﾘﾝｸ(japanese.chosun.com)

680:デフォルトの名無しさん
06/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:デフォルトの名無しさん
06/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:デフォルトの名無しさん
06/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:デフォルトの名無しさん
06/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:デフォルトの名無しさん
06/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:デフォルトの名無しさん
06/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:デフォルトの名無しさん
06/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:デフォルトの名無しさん
06/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
06/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:デフォルトの名無しさん
06/07/19 07:56:54
こんにちは、fortran90での課題が手に負えないので助けてください。
「n×nの実数正方行列Ａの行列式｜Ａ｜を再帰手続きでもとめよ」というもの
ですが、再帰手続きの使い方、使いどころが見当もつきません。３×３までは
できたんですが再帰手続きを使っておらず、それ以上に拡張できません。
長いこと詰まっています。どうか教えてください。

690:デフォルトの名無しさん
06/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:デフォルトの名無しさん
06/07/19 10:48:05
690さん、ありがとうございます。アドバイスをもとに考えてみます。
fortranをはじめたばかりの初心者なので、不躾なんですがよろしければ
なにか具体的な例を示していただければ幸いです。行列を読み込んで添字に
着目して、部分行列を作ってそれから再帰手続きとしようとしているんですが、
再帰手続きの手順がこんがらがりよく理解できません。どうかお願いします。

692:デフォルトの名無しさん
06/07/19 11:19:36
>>689
690 ではない者ですが、余因子展開すると n×n 行列の行列式を
n-1×n-1 行列の行列式で表せると思います。

693:デフォルトの名無しさん
06/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:デフォルトの名無しさん
06/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:デフォルトの名無しさん
06/07/19 15:45:01
＞＞６９３の式では、小行列式を作るのに、あらわにワーク行列を作っているが、
ＦＯＲＴＲＡＮ９０なら、マスキングとかを作って、ワークをあらわに経由しないで、
小行列を引数として再帰関数を呼び出せるような気がする。

手元に紙マニュアルがないので、調べる気がしない。だれか頼む。

でもどっちにしろmask行列を作らないと駄目かな？
意味ないかもｗ

696:デフォルトの名無しさん
06/07/19 15:51:12
687さん、早速お答え頂きどうもありがとうございます！
全体と個々のstepとでloopを回すのでは全く違いますね。
ありがとうございました。

697:デフォルトの名無しさん
06/07/19 16:24:45
692さん、693さん、694さん、695さん、具体例およびアドバイスありがとうございます。
本当に助かります。早速これらのプログラムを参考に自分なりに考えてみます。

698:デフォルトの名無しさん
06/07/20 00:14:15
すいません、fortran77で、
「sin(x)をマクローリン展開し、第ｎ項まで計算したものを、
ｘを0から180度まで30度ずつ変えて表示するプログラムを作成せよ。
但しsin(x)の計算はサブルーチンを使うこと」
という課題が出されたのですが、うまく作れません。
どうかヒントだけでも教えてもらえないでしょうか。

699:デフォルトの名無しさん
06/07/20 03:03:56
>>698
多分マクローリン展開が間違ってるんじゃないかな？

700:デフォルトの名無しさん
06/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
06/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:デフォルトの名無しさん
06/07/20 12:58:47
>>700
SIN(x)=x-x^3/3!+x^5/5!-x^7/7!
これは偶数次の展開が0になって消えているから7次(または8次)の項まで
展開したという理解で良かったっけ?

703:700
06/07/20 13:50:19
>>702
ういうい。
sin は奇関数だから、偶べきの項は出てこない。ゆえに、偶数次の項は消える。

704:困ってます
06/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:デフォルトの名無しさん
06/07/20 16:15:40
>>704
こんなの平均値を求めるのと変わらんぞｗ
本質は１行だｗ
あっちのスレで解いてやったが、今後はこっちで継続せよ。

706:困ってます
06/07/20 16:50:26
あれでやってみたんですが全くできません。。。
なんとか馬鹿でも分かるようにおしえてください！！！！

707:デフォルトの名無しさん
06/07/20 16:57:12
>>706
配列ｆの中に、答えが入っている。
これを出力してグラフを書けばいい。

　＿
/　　\
　　　　\＿/
⇒ｘ軸

　\＿/
⇒ｙ軸　　ｘ＜０

　＿
/　　\
⇒ｙ軸　　ｘ＞０

こんな感じになっている。

708:困ってます
06/07/20 17:10:13
ありがとうございます！！！
でもそれでもいまいち分かってない自分がいます。。。
プログラムをそのままコンパイル出来る形で教えていただけないでしょうか？？？？

709:困ってます
06/07/20 17:13:06
出力したいんですがやり方がやっぱわかないです！！！！！

710:デフォルトの名無しさん
06/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

711:デフォルトの名無しさん
06/07/21 00:18:35
教えてください！！

program oppai
implicit none integer:: saikou,saitei,mark,n
n=0
saikou=0
saitei=100
do
read *,mark
if (mark<0) exit
n=n+1
if (mark<=saitei) then
saitei=mark
end if
if (mark>=saikou) then
saikou=mark
end if
end do
print *,'最高点',saikou
print *,'最低点',saitei
pause
end

このプログラムは、生徒のテストの点数を次々に読み込み、
負の数が入力されたら入力を終了し、その最高点、最低点を出力する
プログラムなんですが、100点満点でないテストでは、うまく動きません。
テストが100点満点かどうかわかっていない場合でも、動くように改良してほしいの
です。。。

お願いします～

712:デフォルトの名無しさん
06/07/21 00:22:43
自然対数の底 e を100桁計算するプログラムを教えてください、
e = 1/0!+1/1! + 1/2! +1/3!+… となります
適当なMを決めて
e = 1/0!+1/1! + 1/2! +1/3!+…+1/M! を計算してやればいいみたいです

ヒントがあります。

これを次のように変形する。

e=(.....(1/M+1)/(M-1)+1)/(M-2)+1)....)/2+1)/1+1
M=80で充分
80! > 10^110
(n! ～ (n/e)^n*sqrt(2*pi*n))
おおまかな求める手順 1.e=1でスタート
2.k=M,M-1,...,1について以下を繰り返す
e←e/k+1 とする

100桁求めるにはM=80で十分

100桁の筆算を計算機上でおこなう
e(0) に整数部分
e(1) に小数第１桁
e(i) に小数大 i 桁

と考え、割り算のプログラムを作る
この100桁の数にたいして kで割るということができればよい。

713:デフォルトの名無しさん
06/07/21 00:23:15

1.234 / 4
-> 1 / 4 --> 0 あまり 1
-> (1*10+2)/4 --> 3 あまり 0
-> (0*10+3)/4 --> 0 あまり 3
-> (3*10+4)/4 --> 8 あまり 2

----> .308 あまり 0.002

e(0) に整数部分
e(1) に小数第１桁
e(i) に小数大 i 桁

e/k (kは1～M=100)

e=(e(0).e(1)e(2)....) をk でわる
amari=0
do i=0,100
jo=amari*10+e(i)
e(i)=jo/k
amari=mod(jo,k)
end do

計算終了後、求めたeの値は次のようにして出力する。

print *, e(0),"."
print '(" ",50I1)',(e(i),i=1,100)
pause
end

お願いします！

714:デフォルトの名無しさん
06/07/21 10:04:01
>>711

saitei=100

これがあるせいで、１００点以上にならない。

saitei = 2**30

こうしておけばおｋ！

715:デフォルトの名無しさん
06/07/21 14:26:10
>>712
これだけヒントをもらっているんだから、もう少しがんばれｗ
かったるいがF77で書いてやった。
１０分で作ったから間違ってても知らんｗ

PROGRAM EXPO70
PARAMETER (ne = 100, n = 80)
INTEGER ie(ne)
ie(1) = 1
DO 10 i = 2, ne
ie(i) = 0
10 CONTINUE
DO 20 m = n, 1, -1
k = ie(1)
DO 30 i = 1, ne - 1
ie(i) = k / m
k = 10 * MOD(k, m) + ie(i + 1)
30 CONTINUE
ie(1) = ie(1) + 1
20 CONTINUE
WRITE(*, '(i1, a1, 99i1)') ie(1), '.', (ie(i), i = 2, ne)
STOP
END

URLﾘﾝｸ(www.math.utah.edu)
値はここで比較しろ。最後の桁が違って居るが、これはたぶんおｋ。

716:715
06/07/21 14:39:13
微妙に改良＋ワーク領域を１０１こ取れば１００桁目まで正しくなる。

PROGRAM EXPO70
PARAMETER (ne = 100, n = 80)
INTEGER ie(ne + 1)
DO 10 i = 1, ne + 1
ie(i) = 0
10 CONTINUE
DO 20 m = n, 1, -1
ie(1) = ie(1) + 1
k = ie(1)
DO 30 i = 1, ne
ie(i) = k / m
k = 10 * MOD(k, m) + ie(i + 1)
30 CONTINUE
20 CONTINUE
ie(1) = ie(1) + 1
WRITE(*, '(i1, a1, 99i1)') ie(1), '.', (ie(i), i = 2, ne)
STOP
END

717:デフォルトの名無しさん
06/07/21 15:42:29
>>716
うは、ne=10000 n=3250で１００００桁までもとまりんぐｗ

718:デフォルトの名無しさん
06/07/22 03:43:44
横35、縦35のマスを考える。そのマス目なかにグループAの人が何人か、グループBの人が何人か、いるものとする。
（各マス目の中には、最大1人しか入ることができない。また、誰もいないマス目もある。）
各人は自分の回りに同じグループの人が多くいた方が幸福だと思っている。そこで「幸福度」を次のように定義する。
自分の場所の上下左右、斜め上下左右、合計8箇所を「自分の回り」とする。この中にいる自分と同じグループにいる仲間の数をMとする。
また、「自分の回り」にいる全ての人の数をNとする。この時、幸福度を
・N=0　なら幸福度 0
・N>0 なら幸福度 M/N
と定義する。

さて、各人は自分の幸福度がある値（以後、しきい値と呼ぶ）より小さいと、現状に対する不満から、ランダムに空いているマス目に移動をしようとする。
ただし、移動しない可能性もある。例えば、上図のように、ある人の「自分の回り」のマス目のうち、
上、下、右斜め下、左斜め上の4箇所が空いており、他のマス目には他の人がいるとする。
この時、この人はそれぞれ確率1/5で
・上に移動
・下に移動
・右斜め下に移動
・左斜め上に移動
・移動しない
という行動をとる。35×35のマスの中には多くの人がいるが、この移動は
・まず、一番上のマスの一番左の人が（もしいれば）移動を試みる。（もちろん、この人の幸福度がしきい値以上なら、移動しない）
・次に一番上のマスの左から2番目の人が移動を試みる。以降、順にその右の人が移動を試みる。
・一番上のマスにいる人の移動が終了したら、次は上から二番目のマスの一番左の人が移動を試みる。
・以下同様。

719:デフォルトの名無しさん
06/07/22 03:44:52
という順番で行われる。
（一般的に言って、この移動を繰り返すと、徐々に同じグループの人が集まってくるようになる。）
以上を前提に次のようなシミュレーションを行う事を考える。
1.初期状態
各マス目に確率 1/3でグループAの人を、確率 1/3でグループBの人を置く。
（これは次のように考えればよい。各マス目でさいころをふり、１か２がでたらグループAの人を置き、
3か4が出たらグループBの人をおく。それ以外はだれも置かない。）
2.各個人は自分の幸福度がしきい値以下であれば、前述の要領で移動を行う。
この移動は前述したように
1.一番上のマスの一番左の人から右へ
2.二番上のマスの一番左の人から右へ
3......
という順番で行われる。
3.2.の移動のプロセスを500回繰り返した後、全員の幸福度の平均値を求める。これを「最終平均幸福度」と仮に名付ける。
（直感的に言えば、幸福度の小さい人が移動するので、初期状態より幸福度の平均値は上がっている事が予想される。）
この「最終平均幸福度」は、各人がランダムに移動するので、シミュレーションを行うたびに異なる可能性がある。
そこで、この値をできるだけ正確に求めるために、この1～3のシミュレーションを20回繰り返し、
「最終平均幸福度」の平均値を求めて出力するものとする。
プログラムの大まかな流れは以下のようになる
1しきい値を決める
2以下を20回繰り返す
2.1マスに初期状態を設定する
2.2幸福度に基づいた各人の移動プロセスを500回行う
2.3「最終平均幸福度」、つまり、2.2を500回行った時点での全員の幸福度の平均を求める。
3プロセス2で行った20回のシミュレーションの各回で求めた「最終平均幸福度」の平均値を求めて出力する。
さて、しきい値がそれぞれ0.3,0.4,0.5,0.6,0.7,0.8,0.9の時の「最終平均幸福度」の平均値を求めて出力するプログラムを作成しなさい。

720:デフォルトの名無しさん
06/07/22 03:46:29
ヒント）各人をランダムに移動させる方法の一例を紹介する。
1.確率1/9で上下左右、斜め上下左右、移動せず、を試みる。
2.移動しない場合は、ここで終了。
移動先が空いていたら、移動して終了。
移動先に別の人がいたら、1.にもどる。

・乱数の発生方法
この問題を解くには、コンピュータに乱数を発生させる必要がある。
これは次のようにして行うことができる。以下は０～１の間の乱数を１０個発生させて、出力するプログラムである。
real::x
integer::i
call random_seed()! これはプログラムの最初に１度だけ呼ぶ
do i=1,10
call random_number(x)! xに 0<=x <1の一様乱数が入る
print *,x
end do
end

プログラムの最初に一度だけ random_seed() というサブルーチンを呼び出す。
あとは、random_number サブルーチンを上記のように呼び出せば、そのたびに、xには 0～1の間の乱数が入る

fortran90です。どなたか手伝ってくれませんか？

721:デフォルトの名無しさん
06/07/22 03:46:30
>>718
壁はどう扱うんだ

722:デフォルトの名無しさん
06/07/22 14:45:52
>>718
ぐだぐだ書いてたら200行ぐらいになったので、↓ここから取ってくれ。
URLﾘﾝｸ(kasamatusan.sakura.ne.jp)

マス目の周りは壁（それ以上移動できない）として処理したつもり。
あんまデバッグしてないから、よく見直した方が良いぞｗ

723:デフォルトの名無しさん
06/07/22 20:08:56
>>722
うは、もう完成させたかｗ

漏れもタラタラやってたら２００行弱ぐらいになった。

閾値が０．７ぐらいで割ときれいに相分離する。

しかし、この問題、結構ややこしいと言うか、移動するときに上から更新していくのだが、
2重占有を避けようとする漏れの糞判定アルゴリズムでは、閾値が大きいと０．８超えると
だんだん下のほうに０が集まってきて、上のほうにａ，ｂがたまってゆくｗｗｗ　

もうちょっと真面目に考えたいｗ

724:デフォルトの名無しさん
06/07/22 22:09:54
>>721
暗黙のルールとして、ます目の周りの壁はそれ以上移動できないようにするそうです

>>722
アリガトウゴザイマス。
ですが、自分は動かせませんでした
もう少し、粘ってみます。

>>723
すみませんが、>>723さんもお願いします。。。

725:718
06/07/22 22:13:46
ついかです。

「最終平均幸福度」は
しきい値が0.3→最終平均幸福度0.65～0.7
　　　　　　 0.6→最終平均幸福度0.9前後
　　　　　　 0.9→最終平均幸福度0.55～0.56
だそうです。

みんな幸せでいいじゃん、とか自分は思うのですが、そんなようにいかないみたいですね

726:723
06/07/23 00:49:38
>>724-725
ういうい。

■大体、最終平均幸福度はそんな感じになっている。
・０．５以下だと、みんな我慢してそんなに動かないで、混じり合ってる。
・０．６～０．７くらいだと、AとBがきれいにすみわけして、仲間同士で島を作る。
・０．８を超えると、動き回りすぎてぎゅうぎゅうと身動きできないところまで行って、
　かえってAとBが混じり合って固定化されてしまう。

■漏れが引っかかっているのは、移動のところで、ここに任意性がある。

>>722氏のを参考に見せてもらったが、前の奴を移動させた新しいマップで
次の位置を移動させている。この場合、右に１個動いたとすると、次の回では
さっき動いたばかりの奴をまた動かすことになる。
（この場合微妙な確率の差で右側にA,Bが溜まるような気がする）

漏れはこれが気になったので、古いマップから新しいマップへ動かすようにしようと
したのだが、この場合古いマップだけ見ていると、行き先の２重占有が起きてしまうので
新しいマップも見なければならないが、そうなると更新時間がずれてしまうので、
非対称性が生じて上のほうにA,Bが溜まってしまう。

まぁ問題文をあらためて読み返すと、更新の方式は一意に決まらないような書き方なので
シンプルな>>722方式が好ましい気がする。

■>>722のプログラムが洗練された感じだったので、漏れももうちょっと清書してから
うｐさせてもらうｗ

727:723
06/07/23 03:51:05
threshold = 0.3: average happiness = 0.6468309
threshold = 0.4: average happiness = 0.7579698
threshold = 0.5: average happiness = 0.8603911
threshold = 0.6: average happiness = 0.9136263
threshold = 0.7: average happiness = 0.9200884
threshold = 0.8: average happiness = 0.5539254
threshold = 0.9: average happiness = 0.5387021
Press any key to continue

上に溜まるといったが、左の間違いだった。配列の行と列を逆にｐｒｉｎｔしていた。
ほの非対称性は更新する順番に依存しているようだ。

一応アルゴリズムを変えて、同じ粒子が２回以上動かないようにした。

あまり清書しなかったｗ　
URLﾘﾝｸ(kasamatusan.sakura.ne.jp)

728:デフォルトの名無しさん
06/07/23 09:07:21
課題が溜まってます。
5×5の行列

1 2 3 4 5
2 3 4 5 6
A=( 3 4 5 6 7)
4 5 6 7 8
5 6 7 8 9

と（括弧は5行分）
1
2
x=(3)
4
5
（括弧は五行分）

1：Aとxのベクトル積
2：A^5
を求めよというものです。

c23456
program enshu2
implicit none
real A(5,5),x(5),y,i,j

data A/1,2,3,4,5,2,3,4,5,6,
p 3,4,5,6,7,4,5,6,7,8,5,6,7,8,9/
data x/5,4,3,2,1/
とりあえずここまで作ってみたものの、計算方法がわかりません。
お願いします。

次ページ