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

1 名前：デフォルトの名無しさん [2006/01/24(火) 09:48:23 ] このスレッドは、他のスレッドでは書き込めない超低レベル、  もしくは質問者自身何が何だが分からない質問を勇気を持って書き込むスレッドです。  FORTRAN使いが優しくコメントを返しますが、  お礼はFORTRANの布教と初心者の救済をお願いします。  37 名前：デフォルトの名無しさん [2006/01/28(土) 19:37:37 ] >>35 そもそも、ｂが式に現れていない。問題を写し間違えていないか？ それとも単に、連立方程式を解いて出たａだけを使うという問題のなのか？ 問題の背景とここまでやった自分の考えを述べなさい。 38 名前：35 [2006/01/28(土) 20:00:48 ] >>37 ええと、問題はこれで合っています。 bは使わないで、問題を解いて出てきたaのみで計算するようです。 行列計算ということなのでガウス消去法を用いて解こうとし それぞれの座標に a(1,1)=1 というように値を当てはめていったのですが、その後計算するにあたって fortranのプログラムにおいてガウス消去法がどのようにやるのか分からず、今に至ります。 39 名前：デフォルトの名無しさん mailto:sage [2006/01/28(土) 20:02:47 ] とりあえずフォートランに多相型っていつはいるんですか？ 40 名前：デフォルトの名無しさん [2006/01/28(土) 20:10:27 ] >>39 polymorphismならFORTRAN2003でそれなりに入っている。 41 名前：デフォルトの名無しさん mailto:sage [2006/01/28(土) 20:10:42 ] >>35 まず、両方の曲線の交点が２つあって、その間の面積は解析的に求めることは簡単ですね。 で、その解析的に出た面積はaが分かれば数値的に出てくると。 そのaは二元連立方程式の解の一つであると。 二元連立方程式はプログラミングしないでも、中学レベルの頭ならすぐ解ける。 どこにプログラミングの必要性があるのか疑問。 42 名前：デフォルトの名無しさん [2006/01/28(土) 20:17:28 ] >>38 いまいち問題設定が見えてこないなｗｗｗ　これが問題のすべてか？背景は無いのか？ ガウスの消去法は、例外処理を考えなければ、手でやるとおりに組めｗ 例外処理はめんどくさい。 ２＊２限定なら、高校で習った公式に突っ込んだほうが早い。例外がdet=0で排除できるし。 というか手でといたほうが早いわけだが。 はっきり言って２＊２をガウスの消去法で解く意味が分からんｗ というわけだから張り切って一般のｎ次行列について作れｗ 43 名前：35 [2006/01/28(土) 20:29:59 ] >>41 これをプログラミングで解いてレポートして提出しろという課題なんですよ。 一応実際に計算して解くことは出来ました。 しかしプログラミングのほうが分からず… >>42 ガウスは使わない方がいいのでしょうか。 となると逆行列を用いての計算のほうがやりやすいのでしょうか？ 44 名前：デフォルトの名無しさん mailto:sage [2006/01/28(土) 20:32:57 ] >>43 じゃあさ、面積を求める部分をSimpson法やいろんな数値積分法で組んでしまうのはどう？ 面積を求める部分をサブルーチンか関数にしてaを渡せばそれなりにプログラミングの形になるし。 45 名前：42 [2006/01/28(土) 20:49:16 ] >>43 ガウスの消去法はアルゴリズムが単純なので、小さめの連立方程式を解くための 練習問題としてはよく使われる。行列がsingularな時や誤差が溜まりやすい状況を、 あまり考慮しなければ、比較的楽に作れるだろう。 素朴には、対角要素が全部０で無いなら、１列目からｎ列目まで、順に各列の対角要素以外を ０にするように行ごとひいていけばいい。手でとくのと同じだ。 対角要素が０や０付近だと例外処理をすることになる。このへんがまんどくさい。 この辺は、問題製作者の意図がどの辺にあるのかで、こだわり度が決まってくるのだが、 この問題設定では、いまいち理解しにくい。 ２＊２での一般形なら単純に公式を使うほうが早いが、n*nの一般サブルーチンを求めているなら、 ガウス法なり何らかの方法が必要になる。 後半では数値積分をすることになるわけだが、問題として全く意図が読めないので、 これだけの情報では、どの方法をとるべきか答えがたい。 46 名前：35 [2006/01/28(土) 21:01:23 ] >>44 うーん、すいません、またｔくの初心者ですのでよくわかりません… 助言は非常に嬉しいのですが…すいません。 >>45 今回はどうやら2*2行列のようです。 それでは今回は公式の方を用いて計算することにします。 まぁ、普通の行列のプログラミングもよくわからないのですが… 後半部分はあれだけです。 タダ単に前半で求めた値を代入して面積を求めるだけでして、背景は無いと思います。 47 名前：デフォルトの名無しさん [2006/01/28(土) 22:46:30 ] >>35 まあ、適当にやれｗ a = 1.0, b = 0.0 交点 x = 0.0, x = e - 1 積分　ｆ（ｘ） = ax - x log(a + x) (a = 1.0　積分範囲　[0.0, e-1]) 解析的な値 -e^2 / 4 + e - 5 / 4 台形公式　短精度　分割数 n=1000 で　こんな感じだ。　 誤差は( (e-1)/1000 )^2〜1e-6位なので、まぁもっともな気がする。 numerical integral 0.378981858 analytic value　　　0.378982186 Program Completed Press Enter to Continue. 48 名前：デフォルトの名無しさん [2006/01/28(土) 22:51:17 ] PROGRAM vip IMPLICIT NONE INTEGER, PARAMETER :: n = 1000 REAL :: x0, x1, s x0 = 0.0 x1 = EXP(1.0) - 1.0 s = trapez_integ(n, x0, x1) PRINT *, 'numerical integral', s PRINT *, 'analytical value ', 5.0 / 4.0 + EXP(1.0)**2 / 4.0 - EXP(1.0) STOP CONTAINS REAL FUNCTION trapez_integ(n, x0, x1) INTEGER, INTENT(IN) :: n REAL , INTENT(IN) :: x0, x1 INTEGER :: i REAL :: s, h h = (x1 - x0) / REAL(n) s = 0.0 DO i = 1, n s = s + 0.5 * ( func(x0 + (i - 1) * h) + func(x0 + i * h) ) END DO trapez_integ = s * h RETURN END FUNCTION trapez_integ ! REAL FUNCTION func(x) REAL, INTENT(IN) :: x REAL, PARAMETER :: a = 1.0 func = a * x - x * LOG(a + x) RETURN END FUNCTION func END PROGRAM vip 49 名前：デフォルトの名無しさん [2006/01/28(土) 22:56:38 ] 訂正 解析的な値 -e^2 / 4 + e - 5 / 4 ↓ e^2 / 4 - e + 5 / 4 正負逆だったｗ 交点が　0.0 と　e-1　なのは　func(0.0) と func(e-1) を書かせて確かめるのもよし。 前半は、２＊２の逆行列の公式でいいだろう。 Ax＝ｂ　 ｘ＝（A^-1）ｂ で一応行列の積も入るし。 50 名前：デフォルトの名無しさん [2006/01/29(日) 12:21:03 ] ひょっとして、数値計算プログラム作るなら Cよりフォートランの方が高速コードを生成しますか？ 51 名前：デフォルトの名無しさん [2006/01/29(日) 13:06:57 ] >>50 昔はそうだった。最近はそれほどで差が無い。 コンパイラ自体がフロントエンド部とバックエンド部に分けて作られるようになり、 中間言語→機械語生成のバックエンド部は共通になったから。 しかしF90以降ならどう考えてもFORTRANの方がバグが少なく、 簡潔明快なソースでやりたいことを表現できる。 52 名前：デフォルトの名無しさん [2006/01/29(日) 13:12:43 ] 速度は？ 53 名前：デフォルトの名無しさん mailto:sage [2006/01/29(日) 13:26:03 ] スパコンのレベルになると、事実上FORTRANしか選択肢がなかったりする でもそこらのパソコンでやるなら、Cで書こうがFORTRANで書こうが大して変わらんね。 g77とかはほとんど内部でf2cを動かして、C言語に変換してからgccでコンパイルしてるし 54 名前：デフォルトの名無しさん [2006/01/29(日) 13:44:16 ] >>52 コンパイラメーカーが同じなら、速度もたいして変わらない。 要するに、ソースを中間言語に直すところでの最適化の違いだけ。 中間言語から機械語へは共通化している。 ただ、FORTRANの方が制約がきつい分コンパイラーの判断できる最適化は Cより高いとされている。 >>53 スパコンに関して言えば、９０年代はじめごろまではベクトル化コンパイラーはFORTRAN だけで、Cコンパイラーがあっても非ベクトル型だった。しかし、その後Cコンパイラーも ベクトル化コードを吐くようになってきた。 55 名前：デフォルトの名無しさん mailto:sage [2006/01/29(日) 14:03:46 ] FORTRAN９０／９５はコードを書いていて気持ちいい。 柔軟性と制約のバランスが取れていてぴったりくる手袋や靴下のよう。 イラつくのは宣言部が重くなりがちな点くらいか。 みなにお勧めする。 56 名前：デフォルトの名無しさん [2006/01/29(日) 14:09:26 ] じゃあ、特に速度を要求する部分を、書けるならアセンブリで、書けないならフォートランで、 んで普段はC++で書くとかって感じですかね 57 名前：デフォルトの名無しさん mailto:sage [2006/01/29(日) 14:26:58 ] MPIのルーチンはCとFortranしか提供されていないが、微妙にCより。 C++は動的バインディングのせいでCにくらべて１０〜２０％遅くなるという話だ。 次期FORTRANにOOP機能が付け加わったらどうなることか。 58 名前：デフォルトの名無しさん [2006/01/29(日) 19:40:26 ] IMPLICIT REAL*8 (A-H,O-Z) WRITE(*,*) 'あ' N=3 COMPLEX*16 A END でコンパイルしようとして $ g77 test3.f test3.f: In program `MAIN__': test3.f:2: WRITE(*,*) 'あ' 1 test3.f:4: (continued): COMPLEX*16 A 2 Statement at (2) invalid in context established by statement at (1) って言われるんですけど、何がいけないんでしょう？ 教えてください 59 名前：デフォルトの名無しさん mailto:sage [2006/01/29(日) 20:22:25 ] >>58 COMLEX文は宣言文だから、最初の実行文よりも前に置かれねばならない。 このばあい、WRITE文の前。 それで怒られている。 60 名前：デフォルトの名無しさん [2006/01/30(月) 10:41:15 ] 2006年の1月から3月までのカレンダーを作れという課題がでました。 ・・・・・・・・・・・ do 20 t=1,3 if(t.eq.2)then b=28 else b=31 end if ・・・・・・・・・ ・・・・・・・・・ write(6,*)" ",t,"月　　　　　　　　" write(6,*)"　日　月　火　水　木　金　土　　" ・・・・・・・・・・ ・・・・・・・・・・ write(6,100)(a(i),i=6,b) 100 format(7i5) ・・・・・・・・・・ ・・・・・・・・・・ ・・・・の部分がわかりません。よろしくお願いします。 61 名前：デフォルトの名無しさん mailto:sage [2006/01/30(月) 11:08:08 ] いまどきフォートラン使ってる人って、どういう立場だなんだ。 趣味？　仕事？　学校？　環境は？　作ってるアプリの種類は？ ちと興味ある。 62 名前：デフォルトの名無しさん [2006/01/30(月) 11:24:39 ] >>61 失敬な！ FORTRANは現役ですよ！いまやモジュール型言語として完成。次はOOPｗ INTELが直接開発・販売しているコンパイラはC＋＋とFORTRANだけですよん。 63 名前：デフォルトの名無しさん mailto:sage [2006/01/30(月) 12:37:58 ] いやだから何作ってるのか教えてよ。現役なのは存じてるから 64 名前：デフォルトの名無しさん mailto:sage [2006/01/30(月) 15:13:53 ] FORTRANを科学技術計算以外で使う人なんているの？ 65 名前：デフォルトの名無しさん mailto:sage [2006/01/30(月) 15:35:34 ] ということはここに居る人みんな科学者！？ 66 名前：デフォルトの名無しさん mailto:sage [2006/01/30(月) 16:43:06 ] >>61, >>63 俺は固体の量子計算。学校の研究です。 アプリっていうとGUIを想像してしまうけど、そうじゃなくて、 コマンドラインで実行するようなプログラム。 >>65 たぶんそうだと思う。 あとガッコウの宿題ができない人々…。 67 名前：デフォルトの名無しさん mailto:sage [2006/01/30(月) 17:29:24 ] Fortran使って構造計算してる土建屋も科学者か？ 68 名前：デフォルトの名無しさん mailto:sage [2006/01/30(月) 21:27:22 ] FORTRAN　ｍｐ３　ｐｌａｙｅｒ www.pbase.com/gadget_square/fortran ﾜﾛｽ 69 名前：デフォルトの名無しさん mailto:sage [2006/01/30(月) 21:55:33 ] fortran mp3 encoder ttp://members.at.infoseek.co.jp/kitaurawa/cgi-bin/wiki.cgi?page=UZURA 70 名前：デフォルトの名無しさん [2006/01/30(月) 21:56:14 ] ところで、fortranに高階関数はいつ入るんですか？ 71 名前：デフォルトの名無しさん mailto:sage [2006/01/30(月) 22:25:00 ] >>70 高階かんすうってなんですか？おしえてかださい 72 名前：デフォルトの名無しさん mailto:sage [2006/01/31(火) 01:46:45 ] 　　　　　　　　　　　 ,,.　 -‐―-　､ 　　　　　　　　　 ／ _..-‐―- 、　　.ヽ 　　　　　　　　 / ／-‐―-..__..ヽ　_..、ヽ 　　　　　　　　.!三三三´　　　 ｀´　　:　i 　　　　　　　　|.!三三　　 __　　￣ _,. .;. i　暗黙の型宣言を守らない 　　　　　　　　|ii　i!"　　´ ｪｪ｀　l´rｪ､ i..!　　　　　　>>60は地獄へ逝く 　　　　　　　　 !ｉ､ i!:::..... -￣r　　ヽ　　.! 　　　　　　　　　!ｉｉ!.|:::::::　　/ ｀−´ヽ　/ 　 _,,..　===＝"/%;;;::::.　　　'∠ﾆﾆ」　 /、　　　　　　　　　　　　_　-‐‐-､ /´　　　　　　 | .%　:::;;;;:　,,,.,,,:. ⌒: ./´　￣＝=-- 　_,,.. -‐‐'"´　_,..-‐ ´ 　　　　　　　 |　 %　　::::::;;;;;;;;;:.""%.／　　　　　　ｒ"　(<二￣￣> 　　　　　　　.|　　%｡ﾟ+。　　。+% ／　　　　　　　｀ヽiﾉ<二.￣￣> 　　　　　　　|　　　%+*ﾟ+。*ﾟ%／　　　　　　　　　（　<二￣￣> 73 名前：デフォルトの名無しさん [2006/01/31(火) 02:26:28 ] >>60　特別にFORTRAN77で書いてやったです。　感謝しやがれです。 PROGRAM calndr C calender for Jan to Mar 2006 INTEGER ia(31) CHARACTER*30 fmt mon1 = 1 C Jan 1 2006 is Sunday (= 1) DO 10 i = 1, 31 ia(i) = i 10 CONTINUE DO 20 it = 1, 3 IF (it .eq. 2) THEN ib=28 ELSE ib=31 END IF WRITE(fmt, '(a, i2, a, i1, a)') 1 '(', 4 * (mon1 - 1), 'x, ', (8 - mon1), 'i4, / (7i4))' mon1 = MOD(ib + mon1 - 1, 7) + 1 WRITE(6, *) ' ', it, '月　　　　　　　　' WRITE(6, *) '　日　月　火　水　木　金　土　　' WRITE(6, fmt) (ia(i), i = 1, ib) WRITE(6, *) 20 CONTINUE STOP END >>73 74 名前：73 [2006/01/31(火) 02:28:02 ] 実行結果です。 1 月　　　　　　 　日　月　火　水　木　金　土 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 2 月　　　　　　 　日　月　火　水　木　金　土 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 3 月　　　　　　 　日　月　火　水　木　金　土 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 Press any key to continue 75 名前：73 [2006/01/31(火) 02:31:04 ] 画面が崩れたので、専ブラ用のアンカーを書いてやるです。　 >>73,75 FORMATの技巧が高度すぎて、>>60はもちろん>>60に出題した どんくさいチビ人間には理解できないだろです！ 76 名前：73 mailto:sage [2006/01/31(火) 02:31:53 ] 専ブラ用のアンカー間違ったです。　 >>73,74 77 名前：デフォルトの名無しさん [2006/01/31(火) 03:14:40 ] >>71 関数を値として引数に取れる関数だよ。 78 名前：デフォルトの名無しさん [2006/01/31(火) 07:00:31 ] 関数のポインタ？ 79 名前：デフォルトの名無しさん mailto:sage [2006/01/31(火) 09:50:36 ] >>77 単にexternalなら今でも出来るが、関数ポインターならFORTRAN2003で可能に。 80 名前：デフォルトの名無しさん [2006/01/31(火) 16:21:20 ] 今日の宿題はまだか？ｗ 81 名前：デフォルトの名無しさん [2006/01/31(火) 22:44:21 ] ところで型推論っていつFortranにはいるんですか？ 82 名前：デフォルトの名無しさん [2006/01/31(火) 23:40:14 ] >>81 型推論ってなんですか？ 83 名前：デフォルトの名無しさん mailto:sage [2006/02/01(水) 00:18:53 ] すみません、質問なのですが。 フォートランで作った実行ファイルとソースファイルから そのソースでcallから呼び出されてる関数 のライブラリ(モジュール？）の場所を探す方法はありますか？ というのも、プログラムを使っていたＵＮＩＸマシンから別のＵＮＩＸマシンにプログラムを移したところ マルチバイト文字が正しくありません。とエラーが出てしまいます。 なのでソースがあるので別のマシンでコンパイル、リンクしようと思ったのです。 そうするとソースの中にcallで呼び出されたサブルーチン記述されていない名前の関数が あったのでリンクのときにこのモジュールを指定しなければいけないんですよね？ その場所がわからずコンパイルでエラーが出てしまいます。 長文で申し訳ないのですがよろしくお願いします。 84 名前：デフォルトの名無しさん mailto:sage [2006/02/01(水) 00:56:50 ] >>83 エラーメッセージを見ていないので状況がいまいちつかめないが、 MAP出力とかのコンパイラ（リンカ、ローダー）オプションで、 必要な情報が得られるかもしれない。 85 名前：83 mailto:sage [2006/02/01(水) 20:00:59 ] >>84 レスありがとうです。 いろいろ悩みましたが ｺﾏﾝﾄﾞ　nm を使ってライブラリを全部grepに渡して関数名を検索したら出てきました。 新しいマシンに移してコンパイルリンクでき、実行も出来ました。 86 名前：デフォルトの名無しさん [2006/02/01(水) 22:34:43 ] 今日は宿題無いのか？ 87 名前：こんにちは [2006/02/01(水) 22:44:03 ] 日付を2つ入力して，その間に何日あるかを書き出すプログラム(Nannichi)を作れ。実行例はつぎの通りである． 日付を２つ入力して下さい． 何年？ 2004 何月？ 7 何日？ 1 何年？ 2005 何月？ 7 何日？ 1 日数は 365日です． これの答えを教えてください　お願いします 88 名前：デフォルトの名無しさん mailto:sage [2006/02/01(水) 22:53:21 ] >>87 >>5-6をみろ。すでにより一般的なうるう年を考慮したサブルーチンが存在している。 それをちょっといじればよい。 質問があれば、日本語でおｋｗｗｗｗｗ 89 名前：デフォルトの名無しさん mailto:sage [2006/02/02(木) 07:55:46 ] >>82 型を宣言しなくても勝手に正しく推論してくれるしくみだよ。 いつはいるの？ 90 名前：デフォルトの名無しさん mailto:sage [2006/02/02(木) 10:02:57 ] >>89 率直に言ってそれは下らないアイデアだと思う。 書きなぐるスクリプト言語ではよくあるが、なんんというかだらしなさのあらわれだ。 どっちにしろ暗黙の型宣言があるし。 91 名前：デフォルトの名無しさん mailto:sage [2006/02/03(金) 01:39:04 ] ROZEN　MAIDEN ・・・　orz 92 名前：デフォルトの名無しさん mailto:sage [2006/02/03(金) 23:41:27 ] ところで、Fortranにproof-carring-codeはいつ入るの？ 93 名前：デフォルトの名無しさん mailto:sage [2006/02/03(金) 23:49:39 ] >>92 proof-carring-codeってなんでちゅか？ばぶー 94 名前：デフォルトの名無しさん [2006/02/04(土) 00:46:34 ] だれかこのプログラミングを教えてください。お願いします。 ｘ=0.1 , 0.2 , 0.3 , 0.4 , 0.5 ｙ=1.228 , 1.005 , 0.823 , 0.674 , 0.552 近似式　ｙ=a*exp(bx) 上の関数とデータで最小自乗法で近似した式を求めよ。 この問題のプログラムがわかりません・・・。お願いします。 95 名前：デフォルトの名無しさん mailto:sage [2006/02/04(土) 01:35:39 ] >>94 PROGRAM vip IMPLICIT NONE INTEGER, PARAMETER :: n = 5 REAL, PARAMETER :: x0(n) = (/0.1 , 0.2 , 0.3 , 0.4 , 0.5 /), & y0(n) = (/1.228, 1.005, 0.823, 0.674, 0.552/) REAL :: y(n), t(n, 2), z(2, 2), zinv(2, 2), v(2) INTEGER :: i y = LOG(y0) t(:, 1) = 1.0 t(:, 2) = x0 z = MATMUL(TRANSPOSE(t), t) v = MATMUL(TRANSPOSE(t), y) zinv = ainverse(z) v = MATMUL(zinv, v) WRITE(*, *) ' least square fit a=', EXP(v(1)), ' b=', v(2) STOP CONTAINS FUNCTION ainverse(x) IMPLICIT NONE REAL, INTENT(IN) :: x(2, 2) REAL :: ainverse(2, 2) REAL :: determinant determinant = x(1, 1) * x(2, 2) - x(1, 2) * x(2, 1) IF (ABS(determinant) < 1.0e-7) STOP ' singular matrix ' ainverse(1, 1) = x(2, 2) / determinant ainverse(2, 2) = x(1, 1) / determinant ainverse(1, 2) = -x(1, 2) / determinant ainverse(2, 1) = -x(2, 1) / determinant RETURN END FUNCTION ainverse END PROGRAM vip 96 名前：デフォルトの名無しさん mailto:sage [2006/02/04(土) 01:36:42 ] >>95 least square fit a= 1.499271 b= -1.998701 1.227659 1.005252 0.8231379 0.6740159 0.5519093 1.228000 1.005000 0.8230000 0.6740000 0.5520000 97 名前：94 mailto:sage [2006/02/04(土) 03:35:31 ] >>95さん　>>96さん ほんとありがとうございます。 月曜にでも学校でやってみます。 また、わからないことがあれば質問させてください。 ありがとうございました。 98 名前：デフォルトの名無しさん [2006/02/04(土) 18:37:30 ] >>95について。 y=a*exp(bx)を次の式に変形し log(y)=log(a)+bx yとxの各値を上の式に代入すると 次の5本の式となる。 log(y1)=log(a)*1+b*x1 log(y2)=log(a)*1+b*x2 log(y3)=log(a)*1+b*x3 log(y4)=log(a)*1+b*x4 log(y5)=log(a)*1+b*x5 ここで、 log(yi)がi番目の要素である列ベクトルをy、 i行1列を1、i行2列をxiとする行列をt、 1番目の要素をlog(a)、2番目の要素をbとする 列ベクトルをvとすると、上の5式は、 y=tv と表せる。 tが1次独立な列から成るので、 vの最小二乗推定は、次の式 (t't)**(-1)t'y となる定理を>>95は使用していると思います。 (注)t'はtの転置を表す。(t't)**(-1)はt'tの逆行列を表す。） 99 名前：デフォルトの名無しさん mailto:sage [2006/02/04(土) 18:56:37 ] ところでＦｏｒｔｒａｎに配列のバウンダリチェックはいつ入るのでしょうか？ 100 名前：デフォルトの名無しさん [2006/02/04(土) 19:05:46 ] >>99 subscript check は、debug 用にoptionで入るのが普通。1960年代からある。 しかし、遅くなるから実際の実行時に使う奴はいない。 この辺はスピード優先の数値計算の価値観を反映している。 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)の部分に変数をいれればいいのですかね？

---

---

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

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