• 締切済み

FDTDの数値計算について

光波の伝搬の計算としてFDTD法を用いる事を考え、Scilabを用いた1次元のFDTD法をまず行うことを考えています。想定する状況は左端に光源があり、空気(屈折率1)を伝搬するという単純なものです。 そこでプログラムを以下のようにしました。 各パラメータは波長ramuda=500*10^-9,dz=ramuda/20,dt=dz/c,透磁率と誘電率mu=ep=1としています。 Hy(1:M)=0; Ex(1:M+1)=0; For n=1:N, Ex(1)=sin(n*dt*c/ramuda); Hy=Hy-dt./mu.*diff(Ex)/dz; Ex(2:M)=Ex(2:M)-dt./ep.*diff(Hy)/dz; end 宇野先生や橋本先生、小舘先生などの本を読む限り条件等はこれでよいはずなのですが、全く伝搬するような計算ができませんでした。数値計算自体慣れていないため、足りない部分があるかと思いますが全く見当がつきませんでした。どうか何が足りていない、何が間違っているなどありましたらご指摘いただきけないでしょうか。

みんなの回答

  • spring135
  • ベストアンサー率44% (1487/3332)
回答No.1

何をやろうとしているのか、よくわかりませんが、波動の伝搬はsin,cosを用いた時点で解を使っているのと同じであり、わざわざ差分化する必要はありません。基本式 ∂^2f/∂^2t=c^2Δf (Δはラプラシアン) 一次元では ∂^2f/∂^2t=c^2∂^2f/∂^2x これを差分化して [f(i+1)-2f(i)+f(i-1)]/(⊿t)^2=c^2[f(j+1)-2f(j)+f(j-1)]/(⊿x)^2 を用いて左端の変位が伝搬していく様子を数値的に追跡するのが伝搬simulationで通常行うことです。 1次元というのにy,zが出てくるのもよくわかりません。 プログラミング以前に目的と手法をよく確認すべきです。

krclight417
質問者

補足

>何をやろうとしているのか、よくわかりませんが、波動の伝搬はsin,cosを用いた時点で解を使っているのと同じであり、わざわざ差分化する必要はありません。 本来であればガウスパルスを用いることを考えていますが,まずは単純な形としてsin波を考えていました. >1次元というのにy,zが出てくるのもよくわかりません。 FDTD法はマクスウェルの方程式を差分化して電場と磁場を交互に計算していく手法なのでy,zが出てくるのは当然かと。というよりもFDTD法をご存知なのでしょうか?

関連するQ&A

  • FDTD法における入射について

    現在FDTD法によるフォトニック結晶光波回路の解析を行っているのですが、入射が思うようにいきません。 僕が行っている解析はEz、Hx、Hy成分で構成されるTM波で、入射のさせかたは、解析領域におけるx=0の位置でEz成分のみについて、 Ez(x,y) = sin((y-y0)π/W) * sin(2πft) y0:波の中心 W :波の幅 f :入射波の周波数 として、時間的にも位置的にも正弦波の連続波として入射させています。Hx,Hyについてはなにも行っていません。  この入射の計算は、FDTDの電磁界を計算する前に入れており、入射の計算を行ったあとに電界と磁界をそれぞれ計算するようにしています。 しかしこのままでは発散していまってので、FDTDの電界の計算は入射位置を除いた(入射位置がEz(1,j)とするとEz(2,j)から)位置から行っています。  一応このやり方で導波路に波は伝搬するのですが、入射初期の波頭は減衰してなくなってしまいます。しかも遠くに伝搬するにしたがって波はどんどん減衰していき、定常応答を得るまでにかなりの時間を要します。波の振る舞いについてはこれで正しいのでしょうか?  あと、時間と位置的にガウス波の孤立パルスも入射させてみたのですが、すぐに減衰して広がってしまって最初の振幅を保ったまま伝搬しません。なにかおかしい所があるのでしょうか?それともこれで正しいのでしょうか?FDTDに関する論文や本などには、入射についてあまり詳しく書いてないので確かめようがありません。どなたかアドバイスしてくれれば幸いです。    質問にまとまりがなく、わかりにくくて申し訳ありませんが、なにかちょっとしたことでもかまいませんので、ご回答をよろしくお願いします。

  • FDTD 遠方界のプログラム

    FDTD法で遠方界を計算するプログラムを作ったのですがコンパイルはうまくいっても実行結果がちゃんと表示されません。プログラムに問題があるのか、大きな値を与えるとbus errorになってしまいます。どなたか遠方界のプリグラムに詳しい方がいらっしゃったら教えて下さい。よろしくお願いします。 プログラム内容 ubroutine far_field implicit none call far_field_yz_plane call far_field_zx_plane call far_field_xy_plane return end subroutine !=================================================================================== ! y-z 平面に対して !=================================================================================== subroutine far_field_yz_plane use consts use fdtd implicit none integer :: i, j, k,ie real(8) :: plane_yz, fac_yz plane_yz = dy*dz fac_yz = plane_yz/(4.0d0*pi*c*dt) ie = nx-6 i = ie+1 do j = 5, ny-5 do k = 5, nz-5 x = (i-0.5d0-ic)*dx y = (j-jc)*dy z = (k-kc)*dz eys = 0.5d0*(ey(i,j,k)+ey(i,j,k+1)) ezs = 0.5d0*(ez(i,j,k)+ez(i,j+1,k)) hys = 0.25d0*(hy(i,j,k)+hy(i,j+1,k)+hy(i-1,j,k)+hy(i-1,j+1,k)) hzs = 0.25d0*(hz(i,j,k)+hz(i,j,k+1)+hz(i-1,j,k)+hz(i-1,j,k+1)) do l = 1, nang time_e = time-dt time_h = time-dt/2.0d0 timesh = -(dirs(l,1)*x+dirs(l,2)*y+dirs(l,3)*z)/c+rf/c tc_e = time_e+timesh tc_h = time_h+timesh m_e = int(tc_e/dt+0.5d0) m_h = int(tc_h/dt+0.5d0) a_e = (0.5d0+tc_e/dt-m_e)*fac_yz b_e = (0.5d0-tc_e/dt+m_e)*fac_yz ab_e = a_e-b_e a_h = (0.5d0+tc_h/dt-m_h)*fac_yz b_h = (0.5d0-tc_h/dt+m_h)*fac_yz ab_h = a_h-b_h wy(l,m_h-1) = wy(l,m_h-1)-hzs*b_h wz(l,m_h-1) = wz(l,m_h-1)+hys*b_h uy(l,m_e-1) = uy(l,m_e-1)+ezs*b_e uz(l,m_e-1) = uz(l,m_e-1)-eys*b_e wy(l,m_h) = wy(l,m_h)-hzs*ab_h wz(l,m_h) = wz(l,m_h)+hys*ab_h uy(l,m_e) = uy(l,m_e)+ezs*ab_e uz(l,m_e) = uz(l,m_e)-eys*ab_e wy(l,m_h+1) = wy(l,m_h+1)+hzs*a_h wz(l,m_h+1) = wz(l,m_h+1)-hys*a_h uy(l,m_e+1) = uy(l,m_e+1)-ezs*a_e uz(l,m_e+1) = uz(l,m_e+1)+eys*a_e enddo enddo enddo ie = 5 i = ie do j = 5, ny-5 do k = 5, nz-5 x = (i-0.5d0-ic)*dx y = (j-jc)*dy z = (k-kc)*dz eys = 0.5d0*(ey(i,j,k)+ey(i,j,k+1)) ezs = 0.5d0*(ez(i,j,k)+ez(i,j+1,k)) hys = 0.25d0*(hy(i,j,k)+hy(i,j+1,k)+hy(i-1,j,k)+hy(i-1,j+1,k)) hzs = 0.25d0*(hz(i,j,k)+hz(i,j,k+1)+hz(i-1,j,k)+hz(i-1,j,k+1)) do l = 1, nang time_e = time-dt time_h = time-dt/2.0d0 timesh = -(dirs(l,1)*x+dirs(l,2)*y+dirs(l,3)*z)/c+rf/c tc_e = time_e+timesh tc_h = time_h+timesh m_e = int(tc_e/dt+0.5d0) m_h = int(tc_h/dt+0.5d0) a_e = (0.5d0+tc_e/dt-m_e)*fac_yz b_e = (0.5d0-tc_e/dt+m_e)*fac_yz ab_e = a_e-b_e a_h = (0.5d0+tc_h/dt-m_h)*fac_yz b_h = (0.5d0-tc_h/dt+m_h)*fac_yz ab_h =a_h-b_h wy(l,m_h-1) = wy(l,m_h-1)-hzs*b_h wz(l,m_h-1) = wz(l,m_h-1)+hys*b_h uy(l,m_e-1) = uy(l,m_e-1)+ezs*b_e uz(l,m_e-1) = uz(l,m_e-1)-eys*b_e wy(l,m_h) = wy(l,m_h)-hzs*ab_h wz(l,m_h) = wz(l,m_h)+hys*ab_h uy(l,m_e) = uy(l,m_e)+ezs*ab_e uz(l,m_e) = uz(l,m_e)-eys*ab_e wy(l,m_h+1) = wy(l,m_h+1)+hzs*a_h wz(l,m_h+1) = wz(l,m_h+1)-hys*a_h uy(l,m_h+1) = uy(l,m_e+1)-ezs*a_e uz(l,m_e+1) = uz(l,m_e+1)+eys*a_e enddo enddo enddo return end subroutine zx平面、xy平面についても同様に書きました。

  • FDTD法を用いた電磁界解析にOpenMPを導入

    FDTD法を用いた電磁界解析にOpenMPを導入するにはどうしたらよいのでしょうか? ちなみに3次元のプログラムを作成しており、電磁界解析のエンジン部分は下記の様になっています。 for( kt=0; kt<=it ; kt++) { /* ※※※forの括弧開始※※※ */ /*********** 1ステップ前の電磁界[jt-1]の値の定式 ***********/ for( jx=0; jx<=ix; jx++ ) for( jy=0; jy<=iy; jy++ ) for( jz=0; jz<=iz; jz++ ) { exg[jx][jy][jz] = ex[jx][jy][jz]; eyg[jx][jy][jz] = ey[jx][jy][jz]; ezg[jx][jy][jz] = ez[jx][jy][jz]; } /*********** 磁界(hx)の計算 ***********/ for( jx=0; jx<=ix; jx++ ) for( jy=0; jy<=iy-1; jy++ ) for( jz=0; jz<=iz-1; jz++ ) { double dey = (ey[jx][jy][jz+1]-ey[jx][jy][jz])/dz; double dez = (ez[jx][jy+1][jz]-ez[jx][jy][jz])/dy; hx[jx][jy][jz] = hx[jx][jy][jz]+dt/amuy[jx][jy][jz]*(dey-dez); } /*********** 磁界(hy)の計算 ***********/ for( jx=0; jx<=ix-1; jx++ ) for( jy=0; jy<=iy; jy++ ) for( jz=0; jz<=iz-1; jz++ ) { double dez = (ez[jx+1][jy][jz]-ez[jx][jy][jz])/dx; double dex = (ex[jx][jy][jz+1]-ex[jx][jy][jz])/dz; hy[jx][jy][jz] = hy[jx][jy][jz]+dt/amuz[jx][jy][jz]*(dez-dex); } /*********** 磁界(hz)の計算 ***********/ for( jx=0; jx<=ix-1; jx++ ) for( jy=0; jy<=iy-1; jy++ ) for( jz=0; jz<=iz; jz++ ) { double dex = (ex[jx][jy+1][jz]-ex[jx][jy][jz])/dy; double dey = (ey[jx+1][jy][jz]-ey[jx][jy][jz])/dx; hz[jx][jy][jz] = hz[jx][jy][jz]+dt/amuz[jx][jy][jz]*(dex-dey); } /*********** 励振 ***********/ for( jx=1; jx<=ix-1; jx++ ) for( jz=1; jz<=iz-2; jz++ ) { ez[jx][80][jz] = sin( 2.0*pi*(kt%100)/100 )*sin(pi*float(jx)/float(ix)); } /*********** 電界(Ez)の計算 ***********/ for( jx=1; jx<=ix-1; jx++ ) for( jy=1; jy<=iy-1; jy++ ) for( jz=0; jz<=iz-1; jz++ ) { double term1 = (1.0-sigmax[jx][jy][jz]*dt/(2.0*epsirx[jx][jy][jz])) /(1.0+sigmax[jx][jy][jz]*dt/(2.0*epsirx[jx][jy][jz])); double term2 = 1.0/(1.0+sigmax[jx][jy][jz]*dt/(2.0*epsirx[jx][jy][jz])); double dhy = (hy[jx][jy][jz]-hy[jx-1][jy][jz])/dx; double dhx = (hx[jx][jy][jz]-hx[jx][jy-1][jz])/dy; ez[jx][jy][jz] = term1*ez[jx][jy][jz]+dt/epsirx[jx][jy][jz]*term2*(dhy-dhx); } ・ ・ ・ この後電界の計算が同様に続きます。 どなたかこのプログラムにOpenMPを導入するにはどうすればよいか教えていただきたいです。

  • 計算式の全微分について

    はじめまして。全微分の問題で式の扱いに困ってしまったものがありまして、お力貸していただければと思い質問させていただきました>< Z=X^n exp(y^m)を全微分!という問題なのですが、途中まで考えてはみたものの。。 dZ=n X^n-1 dexp(y^m) + *** dX^n んーexp(y^m)を微分するとどうなるのか(***の部分) 表記しづらい計算式で申し訳ないのですが、expの扱いがどうもひっかかってます。 よろしくお願いいたします。

  • 積分の計算

    積分について質問です。高橋陽一郎先生の「微分方程式入門」より、積分の計算が一部わからないところがあったので質問させていただきます。 非常に見にくいので、画像を添付しておきますのでそちらを見ていただいたほうが良いと思います。 まず Φ_m(t,s)=∫[s,t]dt_1 ∫[s,t_1]dt_2∫[s,t_2]dt_3…∫[s,t_m-1] A(t_1)A(t_2)…A(tm)dtm と置きます。A(t)はn次正方行列値関数です。 次の文が分かりません。 「こうおくと、Φ_m(t,s)は、 φ= (∫[t,s]A(r)dr)^m = ∫[t,s]dt_1∫[t,s]dt_2…∫[s,t]dt_m A(t_1)A(t_2)…A(t_m) これの右辺において、t_1,....,t_m を大きい順に並べ替えた後に積分したものの1/m! 倍だとo考えることができる」 とありますが、どうやってm!がでてくるのか分かりません。 具体的な計算方法を示してもらえば幸いです。

  • 偏微分方程式の数値解析

    お忙しいところすみません. 現在,FDTD法による数値解析を行っています. FDTD法は空間及び時間微分に対して中心差分を行っていますがこれを4次のルンゲクッタ法で行うことは出来るのでしょうか?ある参考書に以下のような「連立1階微分方程式」を4次のルンゲ-クッタ法で解くというのがありました. dy/dt=f(t,y,z) dz/dt=g(t,y,z) これを4次のルンゲ-クッタ法で解くと y(n+1)=y(n)+1/6*(k1+2*k2+2*k3+k4) z(n+1)=z(n)+1/6*(j1+2*j2+2*j3+j4) ここで k1=Δt*f(t(n),y(n),z(n)) j1=Δt*g(t(n),y(n),z(n)) k2=Δt*f(t(n)+Δt/2,y(n))+k1/2,z(n)+j1/2) j2=Δt*g(t(n)+Δt/2,y(n))+k1/2,z(n)+j1/2) k3=Δt*f(t(n)+Δt/2,y(n))+k2/2,z(n)+j2/2) j3=Δt*g(t(n)+Δt/2,y(n))+k2/2,z(n)+j2/2) k4=Δt*f(t(n)+Δt,y(n)+k3,z(n)+j3) j4=Δt*g(t(n)+Δt,y(n)+k3,z(n)+j3) これをマクスウェルの方程式に対応させ ε∂E(t,x,y,z)/∂t=f(t,H)=∇×H μ∂H/∂t=g(t,E)=-∇×E として出来ないこともない気がするのですが. なぜこの方法を試みようとした理由は現在,通常のFDTD定式化により数値解析を行っていますが ・プログラムは正しい ・始めの部分は正常に計算出来ている ・途中から凹凸が出始め,それが成長して発散している という事態になり,数値的不安定の可能性が強いと思ったからです.オイラー法やルンゲ-クッタ法ではこのような異常振動は進み幅を小さくすれば止まる(ということが数学的に証明されている)ようです. 何卒回答の方,よろしくお願い致します.

  • 磁性 反磁場の計算

    磁性 反磁場の計算 例えば1テスラの磁場を球状の試料に印加したときに磁化が100emu/g だったとき、試料に加わる反磁場はいくらくらいになりますか? Hd = N・M/μ0 (N:反磁場定数 球の場合1/3, M:磁化, μ0:真空の透磁率) という式が教科書に乗っていましたが、磁性の最大の敵である単位にいつも悩まされてしまいます。 よろしくおねがいします。

  • 2重Σを含んだ計算問題

    1/(n^s)={1/Γ(s)}•{¥integral_{0}^{∞} e^{-nt}•t^{s-1} dt} であることを用いて、 Σ_{m=1}^{∞} Σ_{l=m+1}^{∞} {1/{m^k}}•{1/{l^s}} = {1/Γ(k)}•{1/Γ(s)} ¥integral_{0}^{∞} ¥integral_{0}^{∞} {x^(k-1)}•{y^(s-1)}•{1/{e^{x+y}-1}}•{1/{e^y-1}} dx dy を示せという問題について、計算結果があわずに困っています。 (ただし、kは自然数、sは実部が0より大きい複素数、¥integral_{0}^{∞}は0から∞までの積分、Σ_{m=1}^{∞}は0から∞までの和、Σ_{l=m+1}^{∞}はm+1から∞までの和を表すものとします。また、関数の収束性と積分と無限和の順序交換は概知であるとします。) 左辺を計算したときに、計算途中で、 左辺 = Σ_{m=1}^{∞} {1/{m^k}}•{1/Γ(s)} Σ_{l=m+1}^{∞} ¥integral_{0}^{∞} {e^(-ly)}•{y^(s-1)} dy となり、この先を計算しても目標の式と合わず、何度も見直したのですが、どこがよくないのかがよくわかっていません。(特に、e^(-ly)の部分がうまく捌けません。) もしもわかられる方がおられれば、お教え頂けないでしょうか?

  • リアクトルのギャップ計算

    ACリアクトルを実際に試験したデータをギャップ計算式にあてはめて、比透磁率urを 算出しようとしたのですが、値がマイナスになってしまいます。 L/n^2=(lc/uru0S+Is/u0S)^-1より インダクタンス0.038H巻数198TカットコアSL1700磁路0.451M ギャップ0.0044M 0.038/198^2=(0.451/ur4πx10^-7x0.00284+0.0044/4πx10^-7x0.00284)^-1 ^-1より上記の式を分子分母を逆にして計算していくと値がマイナスになってしまいます。 アドバイスよろしくお願いいたします。

  • fortran FFT

    FDTD法を用いた計算の中で遠方界を計算するのにフーリエ変換が必要になり以下のプログラムを作ったのですがコンパイルの際に以下のエラーが出てしまいます。文字の定義もいろいろと変えてやってはいるのですが、なかなかうまくいきません。どなたかわかる方がいらっしゃったら教えて下さい。よろしくお願いします。 plane_yz = dy*dz fac_yz = plane_yz/(4.0d0*pi*c*dt) ie = nx-6 i = ie+1 do j = 5, ny-5 do k = 5, nz-5 x = (i-0.5d0-ic)*dx y = (j-jc)*dy z = (k-kc)*dz eys = 0.5d0*(ey(i,j,k)+ey(i,j,k+1)) ezs = 0.5d0*(ez(i,j,k)+ez(i,j+1,k)) hys = 0.25d0*(hy(i,j,k)+hy(i,j+1,k)+hy(i-1,j,k)+hy(i-1,j+1,k)) hzs = 0.25d0*(hz(i,j,k)+hz(i,j,k+1)+hz(i-1,j,k)+hz(i-1,j,k+1)) do l = 1, nang time_e = time-dt time_h = time-dt/2.0d0 timesh = -(dirs(l,1)*x+dirs(l,2)*y+dirs(l,3)*z)/c+rf/c tc_e = time_e+timesh tc_h = time_h+timesh m_e = int(tc_e/dt+0.5d0) m_h = int(tc_h/dt+0.5d0) a_e = (0.5d0+tc_e/dt-m_e)*fac b_e = (0.5d0-tc_e/dt+m_e)*fac ab_e = a_e-b_e a_h = (0.5d0+tc_h/dt-m_h)*fac b_h = (0.5d0-tc_h/dt+m_h)*fac ab_h = a_h-b_h wy(l,m_h-1) = wy(l,m_h-1)-hzs*b_h wz(l,m_h-1) = wz(l,m_h-1)+hys*b_h uy(l,m_e-1) = uy(l,m_e-1)+ezs*b_e uz(l,m_e-1) = uz(l,m_e-1)-eys*b_e wy(l,m_h) = wy(l,m_h)-hzs*ab_h wz(l,m_h) = wz(l,m_h)+hys*ab_h uy(l,m_e) = uy(l,m_e)+ezs*ab_e uz(l,m_e) = uz(l,m_e)-eys*ab_e wy(l,m_h+1) = wy(l,m_h+1)+hzs*a_h wz(l,m_h+1) = wz(l,m_h+1)-hys*a_h uy(l,m_e+1) = uy(l,m_e+1)-ezs*a_e uz(l,m_e+1) = uz(l,m_e+1)+eys*a_e enddo enddo enddo コンパイル後timesh = -(dirs(l,1)*x+dirs(l,2)*y+dirs(l,3)*z)/c+rf/c 1 Error: Unexpected array reference at (1) In file far_field.f90:48 wy(l,m_h-1) = wy(l,m_h-1)-hzs*b_h 1 Error: Unexpected array reference at (1) In file far_field.f90:49 wz(l,m_h-1) = wz(l,m_h-1)+hys*b_h 1 Error: Unexpected array reference at (1) In file far_field.f90:50 uy(l,m_e-1) = uy(l,m_e-1)+ezs*b_e 1 Error: Unexpected array reference at (1) In file far_field.f90:51 と続きます。