• 締切済み

【fortran】フーリエ級数について

こんばんは。 早速ですが、質問させて頂きます。 学校の課題で一次元熱伝導方程式のプログラムを組みなさいという課題が出たのですが、入力、出力のプログラムが上手く書けずに困っています。 自分で何回もやってみたのですが、「ループがベクトル化されました」としか出ず、どこが間違っているのかも分からない状況です。 詳しい方、教えて頂ければ幸いです。 下に自分が作成したプログラムを載せておきます。 program heat implicit none integer i, n, nmax real(8) t0, t1, d, time, a, h, p, q, an0, kn, bn, en integer, parameter :: id = 100 real(8), parameter :: pi = 3.1415926535d0 real(8), parameter :: lambda = 2.5d-6, rho = 2500.d0, c = 2.35d-7 real(8), parameter :: expmax = 50.d0 real(8) x(0:id), t(0:id) ! ! Data input part (t0,t1,d,a,time,nmax) ! t0=300.0d0 ! degrees t1=200.0d0 ! degrees d =100.0d0 ! m time=1.0 ! year time = time *365.0d0*24.0d0 ! year --> hour nmax=1000 a = 1.0d6 ! Area of rock surface [m^2] h=d/id do i=0,id x(i)=i*h t(i)=t1 end do p=0.0d0 q=0.0d0 do n=1,nmax,2 an0=4.0d0/(n*pi)*(t0-t1) kn=n*pi/d bn=lambda*kn*kn/(rho*c) if(bn*time.le.expmax) then en=exp(-bn*time) else en=0.0d0 endif do i=0,id t(i)=t(1)+an0*sin(kn*x(i))*en enddo p=p+an0*kn*en q=q+an0*kn/bn*(1.0d0-en) enddo p=lambda*a*p q=lambda*a*q open(20,file=‘output.dat’) write(20,*) ‘x(i),t(i),p,q' = ‘, ‘x(i),t(i),p,q' close(20) stop end program heat

みんなの回答

回答No.2

>write(20,*) ‘x(i),t(i),p,q' = ‘, ‘x(i),t(i),p,q' i の値は何? 上の方で i のループが終わっているから、配列の範囲外では?

  • Tacosan
  • ベストアンサー率23% (3656/15482)
回答No.1

Fortran で open(20,file=‘output.dat’) とか書けたっけ? 引用符, これでいいの? 「上手く書けない」ってのが, 何がどう「上手く書けない」のかわからんけど....

関連するQ&A

  • フーリエ級数の展開

    An=2/4∫(0→4)i(t)sin(nωt)dt Bn=2/4∫(0→4)i(t)cos(nωt)dt i(t)=0:0≦t≦1 =(2-t)Im:1≦t≦2 =0:2≦t≦3 =(t-4)Im:3≦t≦4 以上のような関数を展開したとき、 An={ {1/(nω)}sin(2nω)-{2/(n^2ω^2)}sin(5nω/2)sin(nω/2)}Imsin(nω) Bn={ {1/(nω)}cos(2nω)-{2/(n^2ω^2)}sin(nω/2)cos(5nω/2)}Imsin(nω) となるそうなのですが、いくら計算してもこのような形になりません。 三角関数同士の積の形にすら辿り着けず苦戦しております。 計算過程や注意点を教えていただければ幸いです。

  • フーリエ級数解析によって位相差から時間差を求めたい。

    フーリエ級数解析によって位相差から時間差を求めたい。 異なる時間に発生する2つの信号S1,S2に対して,各々,フーリエ級数解析を行い,求められる各々の三角級数の位相差を時間差として求めたいと考えています。 フーリエ係数をAn,Bnとすると,各々の位相P1,P2は P1=atan(Bn1/An1) P2=atan(Bn2/An2) で求められると思います。 よって,位相差は, |P1-P2| で求められ,この位相差は角速度ラジアンで与えれるために時間(秒)に変換する必要があると解釈しています。 (今回は2πを3.2秒とする。) したがって,時間差は, time(秒)=(3.2*|P1-P2|)/2π で求められる。。。 と解釈しているのですが,理論値とは全く異なる値が求まります。 何が,いけないのかがわからずにおります。 詳細,ご教授いただければ,嬉しく思います。 よろしくお願いいたします。

  • 矩形波『Dutyアリ』のフーリエ級数で周波数は?

    矩形波『Dutyアリ』のフーリエ級数展開が >矩形波f(t)の1周期を[-π,π]、duty=d (0<d<1)とし、f(t)を >f(t)=1, |t|≦dπ >f(t)=0, dπ<d≦π >で与えるとすると、f(t)は偶関数となるから >フーリエ係数a0,an,bnは >bn=0 (n=1,2,3,...) >a0=(1/π)∫[-dπ,dπ] 1 dt=2dπ/π=2d >an=(1/π)∫[-dπ,dπ] 1*cos(nt)=2sin(nπd)/(nπ) (n=1,2,3,...) > >∴f(t)=d+Σ[n=1,∞] {2/(nπ)}sin(nπd)cos(nt) >=d+(2/π)[sin(πd)cos(t)+(1/2)sin(2πd)cos(2t) >+(1/3)sin(3πd)cos(3t)+(1/4)sin(4πd)cos(4t)+ … >+(1/n)sin(nπd)cos(nt)+ …] >と展開できます。 までは理解できるのですが, これに 周波数[Hz]を組み込むには どのような式になるのでしょうか?

  • 質問(Fortran)について

    おはようございます。 下記のFortranでここから~ここまでと---で囲んでいる 部分の意味が分かりません。 乱数Pを発生させてPがEE以下になったら、140へ行けとあるのですが。。 それと(1)の部分をCに書き換えるとどうなるか、知っている方教えてください。 よろしくお願いいたします。 C *---------------------------------------------------------* C サブルーチンWGEN C *---------------------------------------------------------* SUBROUTINE WGEN(EM,R,NN,IR,ACC,ND,DT,AMAX,VMAX,MXCYCL,ERR,UW1, * UW2) C implicit real*8(a-h,o-z) COMPLEX*16 C(4096) DIMENSION ACC(ND),UW1(ND),UW2(ND) DIMENSION E(33),X(33),EE(33) DIMENSION PDIF(2046),PHI(2049),F(2049),T(2049),SV(2049),H(1), * RES(2049,1),RR(2049) PARAMETER (PI2=6.283185) DATA DX/0.03125/,H0/0./,H/0.05/,IR/101/ C C ------------------- ここから ------------------- DO 150 K=1,NN2-2 P=RAND2(IR) ← (1) C DO 130 J=2,33 IF(P.LE.EE(J)) GO TO 140 C 130 CONTINUE C 140 PDIF(K)=-(X(J-1)+(P-EE(J-1))/(EE(J)-EE(J-1))*DX)*PI2 C 150 CONTINUE C ------------------- ここまで ------------------- C *---------------------------------------------------------* C ファンクションRAND01 C *---------------------------------------------------------* REAL*8 FUNCTION RAND2(I) C INTEGER*4 L,C,T30 REAL MU PARAMETER(L=843314861,C=453816693,T30=2**30,MU=2.0**31) C I=L*I+C IF(I.LT.0) I=(I+T30)+T30 RAND2=REAL(I)/MU END C

  • フーリエ変換です. よろしくお願いいたします.

    フーリエ変換です. 周期T=2πのf(t)=1-t/π (0≦t≦2π) があるとき, tn=nT/N (0≦n≦N-1), ωp=2πp/T (0≦p≦N-1)を用いて f(tn)=Σ[p=0~N-1] {X(ωp)e^(iωp・tn)} X(ωp)=1/NΣ[q=0~N-1]{f(tq)e^(-iωp・tq)} このとき, X(ωp)を求めよ. という問題で, X(ωp)=1/NΣ[q=0~N-1]{f(tq)e^(-iωp・tq)} から X(ωp)=1/NΣ[q=0~N-1]{(1-q/π)e^(-iωp・tq)} =1/NΣ[q=0~N-1]{(1-q/π)e^(-i2πpq/N)} =∫[q=0~N-1]{(1-q/π)e^(-i2πpq/N)} と式変形をしてX(ωp)を出そうとしたのですが,うまくいきませんでした. 解き方が間違っているのでしょうか. よろしくお願いいたいします.

  • FortranをCに書き直すにはどうしたらいいでしょうか?

    Foratranで書かれている下記のコードをCのコードに置き換えたいのですが、分かりません。 お願いします=(^v^)= C *---------------------------------------------------------* C サブルーチンWGEN C *---------------------------------------------------------* SUBROUTINE WGEN(EM,R,NN,IR,ACC,ND,DT,AMAX,VMAX,MXCYCL,ERR,UW1, * UW2) C implicit real*8(a-h,o-z) COMPLEX*16 C(4096) DIMENSION ACC(ND),UW1(ND),UW2(ND) DIMENSION E(33),X(33),EE(33) DIMENSION PDIF(2046),PHI(2049),F(2049),T(2049),SV(2049),H(1), * RES(2049,1),RR(2049) PARAMETER (PI2=6.283185) DATA DX/0.03125/,H0/0./,H/0.05/ C C DO 150 K=1,NN2-2 P=RAND2(IR) C DO 130 J=2,33 IF(P.LE.EE(J)) GO TO 140 C 130 CONTINUE C 140 PDIF(K)=-(X(J-1)+(P-EE(J-1))/(EE(J)-EE(J-1))*DX)*PI2 C 150 CONTINUE C C *---------------------------------------------------------* C ファンクションRAND01 C *---------------------------------------------------------* REAL*8 FUNCTION RAND2(I) C INTEGER*4 L,C,T30 REAL MU PARAMETER(L=843314861,C=453816693,T30=2**30,MU=2.0**31) C I=L*I+C IF(I.LT.0) I=(I+T30)+T30 RAND2=REAL(I)/MU END C

  • fortran77教えてください

    fortran77のプログラムについての質問です。 次のプログラムを実行するとどのような結果になるか教えてください REAL A,B,C,D,E,F A=7.0 B=5.0 CALL WASA(A,B,C,D) CALL WASA(C,D,E,F) WRITE(*,*)E,F STOP END SUBROUTINE WASA(P,Q,R,S) REALP,Q,R,S R=P+Q S=P-Q RETURN END

  • 漸化式と極限の問題

    p,qは正の有理数で、√qは無理数であるとする。自然数nに対し、有理数An,Bnを(p+q)n乗=An+Bn√qによって定める。 (1)(p-√q)n乗=An-Bn√qを示せ。 (2)An   ― ――→ √q を示せ。    Bn (n→∞) 出典:03大阪市大 この問題の解答を教えていただきたいです。 (1)は数学的帰納法を使えば解けるのは分かっているのですが、n=1がどう説明すれば成立を示せるのかが分からず、困っています。 よろしくお願いします。

  • 等差数列の共通項

    等差数列{An}と{Bn}があります。この2つの等差数列の共通項を並べてできる新しい数列の一般項を求める問題です。 {An},{Bn},それぞれ一般項が, An=8n-2, Bn=6n+2 です。 また,それぞれの項を少し書き出すと, {An}:6,14,… {Bn}:8,14,… と,共通項の最小値が14であることが分かります。 ここで,{An}の第p項と{Bn}の第q項が等しいとすると, Ap=Bq であるので, 8p-2=6q+2 となります。 よって, 4p=3q+2 となり,変形して, 4(p-2)=3(q-2) と表されます。 ここまではよいのですが,次のkの置き方について,問題集の回答を見たのですが,いまいちよく分かりません。以下はその解答です。 「4と3は互いに素であるので,kを自然数として,  p-2=3(k-1), q-2=4(k-1)」 何故ここで(k-1)なのでしょうか?kではいけない理由は何でしょうか? どなたか分かる方,教えてください。

  • 「二つのセルの条件によってセルの色替」をマクロで実施するには?

    エクセルファイルの任意セルのコメントによって、あるセルを色付けする方法をご教示いただけないでしょうか。 列 An Bn Cn Dn  En Fn Gn Hn  In Jn Kn Ln  (n=1,2,3,・・・n) An~Hnのセルには「使用」または「未使用」のコメントが入力されています。 InはAn,En、JnはBn,Fn、KnはCn,Gn、LnはDn,Hnのコメントによって色付けします。 <やりたいこと> 1. In~Hnが空欄なら、2,3,4の処理は行わない。 2. En=「使用」かつAn=「未使用」 → In:黄で色付け   Jn,Kn,Lnはそれぞれの対象セルの状態で色付け 3. En=「未使用」かつAn=「使用」 → In:シアンで色付け   Jn,Kn,Lnはそれぞれの対象セルの状態で色付け 4. En=An=「使用」 → In:ピンクで色付け   Jn,Kn,Lnはそれぞれの対象セルの状態で色付け  *行数nは、シートによって様々です。最大512行です。 とりあえず条件付書式を設定して対応したのですが、マクロを使用すれば簡便化できると聞きました。しかし初心者の為、マクロ知識が乏しく判りません。同じ処理を施したいエクセルファイルが多数あり、なんとか効率よく作業ができないかと思案しています。どなたかマクロ処理をご教示していただけないでしょうか。宜しくお願い致します。

専門家に質問してみよう