- ベストアンサー
Doループとファンクション（fortran)
課題で初めてプログラムに取り組んでいるものです。よろしくお願いします。 台形則を使って楕円の面積（の４分の１）を求めるプログラムを作ったのですが、ファンクションで関数を指定するとうまくいきません。 具体的には func(x)=4.0*sqrt(1.0-x**2.0) という関数で、0<=x<=1を１０００万分割して面積を近似するという課題です。 本来ならfunc(x)は徐々に減少し、x=1で０になるはずが、実際に走らせてみると、それより手前で０になってあとはNaNがでます。 func(x)以外（x座標の分割の仕方等）はうまく行っており（確認済み）、function形式をやめてメインのプログラムに直接式を書き込んだらうまくいきました。 自分ではなんでこんなことが起こるのか皆目見当が付きません。ご教授お願いしますm(_ _)m
- silk-bombaye
- お礼率100% (2/2)
- その他（プログラミング・開発）
- 回答数2
- ありがとう数4
- みんなの回答 （2）
- 専門家の回答
質問者が選んだベストアンサー
a，bそれぞれを不動小数点数とすると、 a^b = e^(b*log(a)) として計算されます。底はeです。 この計算だけでも、EXP(X)の計算，ALOG(X)の計算，乗算を経由するたびに誤差が積もリます。そのうえ、単精度浮動小数点数と倍精度浮動小数点数の精度の違いで、誤差がいっそう積もりますので、べき計算は"X**2"とすべきでしょう。 関数で計算する場合は、ご自身で倍精度計算に変更してみるのはいかがでしょう。 4.0*sqrt(1.0-x**2.0) ⇒DBLE(4.0)*DSQRT(DBLE(1.0)-X**2) 引数Xは倍精度とします。関数FUNNCも倍精度宣言すれば、メイン側で計算するのと同じ結果が得られるでしょう。
その他の回答 (1)
- ultraCS
- ベストアンサー率44% (3956/8947)
課題だと答えを出すのは為にならないので たぶん、誤差の問題ですね。 何も宣言しないと、XやFUNC、1.0、2.0は単精度になっているはずです。単精度だと、有効数字の範囲が、1000万をフルに表現できるだけありません。 ではどうするかは、ちょっと考えてみてください。 数値計算に取り組むなら、まず、誤差について勉強してみましょう。 なお、余談ですが、べき乗は処理系にもよりますが、X**2.0とX**2、X*Xでは内部の処理が違います、べき数を浮動小数点にすると、一度対数を呼ぶので計算課程が複雑になり、精度、速度が下がることがあります。2乗であれば、x**2かX*Xと書きましょう。
お礼
有難うございました!! 数値を倍精度にする方法はまだ知らないので、自分で探してみます(^^)
関連するQ&A
- fortran90 引数で渡された関数の呼び出し
fortran90を始めて間もない者です。 メインプログラムより呼び出しているサブルーチンにユーザ関数を渡しています。 このサブルーチンを自前で作成するのが目的です。 引数で受け取った関数を、自前の別関数より呼び出すにはどうしたら良いでしょうか。 Cであれば関数のポインタをグローバルな変数にセットしてやれば可能だと思いますが・・・ !-------------------------------- subroutine sub(func1, a) real::a interface real function func1(x) real::x end function func1 real function func2(x) real x end function func2 end interface call sub2(func2, a) write(*,*) a return end subroutine sub !-------------------------------- real function func2(x) real::func2, x ! ここでfunc1を呼び出したい ! func2=func1(x) end function func2 !-------------------------------- subroutine sub2(funca, a) real::a interface real function funca(x) real::x end function funca end interface a=funca(10.) return end !-------------------------------- program main external func real a call sub(func, a) write(*,*) a end program !-------------------------------- function func(x) real func, x func=2.*x*x end function func
- 締切済み
- その他（プログラミング・開発）
- fortranのプログラムについて
fortranでプログラムを作成する課題が出たのですが、やり方が全くわかりません。 何を使ってどういう順番で組めばよいか教えてください。 課題は以下のとおりです。 x^2+y^2<1を満たす領域の面積を、以下の方法で近似計算するプログラムを作成せよ。 はじめは、領域を0<x<1、0<y<1と設定し、x、y方向にそれぞれn等分に分割し、多数の正方形を用意する。 分割されたそれぞれの正方形の面積は、みな等しく1/n^2である。 各正方形の重心位置（xi,yj）は、 xi=1/2n+i/n i=0～n-1 yj=1/2n+j/n j=0～n-1 であらわされる。 個々の正方形について、その重心位置（xi,yj）が求めたい領域に含まれるかどうかを判断し、この総数から近似面積を求める。最後に、この値を４倍したもの（つまりx^2+y^2<1の近似面積）を出力することとする。 なお、i,j,nはinteger（8）で、x,y,求めたい面積等はreal（8）を用い、n=100000（nはread文で入力する形でよい）で計算することとする。 よろしくお願いします。
- 締切済み
- その他（プログラミング・開発）
- fortranでのプログラムについて
初歩的な質問失礼します func = x**3-3*(x**2)-2*x+14 に対して real function ymax( minx, maxx) real, intent(in) :: minx,maxx real :: mx1,mx2 mx1 = (6.0+sqrt((-6.0)*(-6.0) - 4.0*3.0*(-2)))/(3.0*2.0) mx2 = (6.0-sqrt((-6.0)*(-6.0) - 4.0*3.0*(-2)))/(3.0*2.0) ! mx2は極大 mx1は極小 if (minx <= mx2 .and. mx2<maxx ) then if (func(maxx)<func(mx2) .and. func(minx)<func(mx2)) then ymax=func(mx2) elseif (func(mx2)<func(maxx) .and. func(maxx)>func(minx)) then ymax=func(maxx) else ymaxx = func(minx) endif else if (func(minx)>func(maxx)) then ymax=func(minx) else yamx=func(maxx) endif endif end function ymax こういう感じで関数の最小最大を決めてグラフを画面上に変換して表示しようとするのですがどうも最大値がうまく計算されません。 どこに問題があるのでしょうか？
- ベストアンサー
- その他（プログラミング・開発）
- FORTRANをご存じの方がいらっしゃいましたら、間違いを指摘していただきたいです。
大学で現在FORTRANのシンプソン公式の課題が出ています。 自分で考えて下記プログラムを作ったのですが、どうしても 正しい答えが出力されません。 正しい答えとしては、分割数４とか５とかで πに限りなく近い値になるそうです。 （下に添付した物は、分割数を4にしています。） もしFORTRANをご存じの方がいらっしゃいましたら どこがどのように間違っているのか、ご指摘いただけないでしょうか？ よろしくお願いいたします。 問題 4/(1＋x＾2)を、積分範囲：0～1 を、シンプソン公式を使って求めよ。 区間分割数については、各自適当に与え 分割数が多くなると精度が良くなる事を確認せよ。 シンプソン公式を使うと台形公式を使う場合より 少ない分割数で解（＝π）に近づくか確認せよ。 ＊＊＊ここから＊＊＊ c numerical value integral c trapezoid formula c *** main program *** real a,b parameter (n=4) a=0 b=1 call simpson(a,b,n,sumf) write(*,*) 'sumf=',sumf end c *** end of main program *** c *** sub program *** subroutine simpson(a,b,n,sumf) dimension f(1000) m=2*n h=(b-a)/float(m) do 10 i=1,m+1 x=a+h*float(i-1) f(i)=func(x) 10 continue sumf=h/3.0*(f(1)+f(m+1)+4*f(m)) do 20 i=2,m sumf=sumf+4.0*(h/3)*f(2*i+1) 20 continue do 30 i=2,m sumf=sumf+2.0*(h/3)*f(2*i) 30 continue end c *** end of sub program *** c *** function sub program *** function func(x) func=4.0/(1.0+x**2) return end c *** end of function sub program *** ＊＊＊ここまで＊＊＊ もしお分かりになる方がいらっしゃいましたら、 お手数ですがご指導いただけないでしょうか？ お願いいたします。
- ベストアンサー
- その他（プログラミング・開発）
- fortran sqrtコンパイルエラー
プログラムを書き直していて以下の様なエラーが出てしまいます。 以前sqrt関数を書いたときは以下の様なエラーは出てこなかったのですが、最近やたらFunctionのエラーがよく出ます。コンパイルに問題があるのでしょうか？自分のコンパイルを調べるにはどうしたらわかるのでしょうか？よろしくお願いします。 In file init.f90:52 rf = sqrt((nx-5)**2+(ny-5)**2+(nz-5)**2) 1 Error: Type of argument 'x' in call to 'sqrt' at (1) should be REAL(4), not INTEGER(4) In file init.f90:52 rf = sqrt((nx-5)**2+(ny-5)**2+(nz-5)**2) 1 Error: Function 'sqrt' at (1) has no implicit type
- ベストアンサー
- その他（プログラミング・開発）
- fortran90での三重積分
fortran90での台形公式を用いた三重積分について悩んでいます。 台形公式を用いた定積分は色んなサイトを見て何となく理解しました。 三重積分は恐らく、３つのfor文の入れ子によるものだと予想できますがソースコードがなかなか記述できません。 被積分関数（例えばf(x)=x^3+y+z^2）をfunctionで定義したいのですが、どなたか御教授願います。
- 締切済み
- その他（プログラミング・開発）
- fortran77
プログラムを作っているのですが、 implicit noneを付けた場合の 関数ｆ（ｘ、ｙ）の宣言の仕方がわかりません。 ググッていくつか出てきたのを試してみたのですが、error文が出てきて、 実行してもうまくいきません。 real x,y,f(x,y) real x,y function f(x,y) real function f(x,y) などです。 implicit noneを付けないとうまくいくのですが・・・・・・ 今後のためによろしくお願いします。
- ベストアンサー
- その他（プログラミング・開発）
- fortranのプログラム
fortranのプログラム 現在、fortranの勉強をしております。 そこで、質問があるのですが、 ある関数ｆ(x,y,z)の座標（x,y,z）の値がデータとして与えられているとき、 S=10+f をfortranで計算したいと考えております。 ただ、関数ｆは複数（ｆ１、ｆ２、ｆ３）あり、次々とｆに代入してＳを計算したいのですが、どのようにプログラムしたらいいか思いつきません。 どなたか、ヒントだけでもいいので、教えてください。 ちなみに、私が考えたプログラムは(下のプログラムはポイントだけ書いてあります。endやその他関係ないと思われるところは省いております。) do 100 k=1,3 S=S+fk(x,y,z) continue function f1 f1(s,t,u)=・・・ return end f2(s,t,u)=・・・ return end 使用しているバージョンは、fortran77(本当は90を使っているのですが、77だけで書いています)です。
- ベストアンサー
- その他（プログラミング・開発）
お礼
整数以外の結果が出る計算をするときは全ての数を浮動小数点にしなければならないと思っていました(^^;; 誤差というものをもっと学ぶ必要がありますね・・・。 数値の倍精度指定の仕方まで教えてくださって有難うございました!!