- ベストアンサー
差分法
単振動の運動方程式 x"+x＝0 を差分法で解こうと思っています。 おそらく差分方程式は、 (x[n+1]-2x[n]+u[n-1])/dt^2=u[n] で良いと思われるのですが（違ってたら指摘してください）、初期条件 x(t=0)=aとx'(t=0)=b の与え方が良く分かりません。 また、sin(ωt)の強制振動を与えたときの差分方程式はどのようになるのでしょうか？ よろしくお願いします。
- alha
- お礼率50% (22/44)
- その他（学問・教育）
- 回答数5
- ありがとう数1
- みんなの回答 （5）
- 専門家の回答
質問者が選んだベストアンサー
そもそもこの方程式の解はは振動するものであってもう少しで発散するものです 近似式が発散しないようにするためにはｈをかなり小さくしなければなりません そのとき発散のスピードは（１＋ｈ＾２）＾（ｎ／２）です この値が１に比べてあまり大きくないとみなすことができるｎの範囲が近似の結果が意味をなすｎの範囲です より大きなｎにおける近似結果を利用する場合にはそれに応じてｈを小さくしなければなりません 元の式が安定と不安定の間になくて安定であればこのような不都合はありません 一応更正版を出しますが勿論式は保証できません （ｘ［ｎ＋２］－２・ｘ［ｎ＋１］＋ｘ［ｎ］）／ｈ＾２＋ｘ［ｎ］＝ｕ［ｎ］ すなわち ｘ［ｎ＋２］＝２・ｘ［ｎ＋１］－（ｈ＾２＋１）・ｘ［ｎ］＋ｈ＾２・ｕ［ｎ］・・・（＊） を初期条件 ｘ［０］＝ａ と （ｘ［１］－ｘ［０］）／ｈ＝ｂすなわちｘ［１］＝ｂ・ｈ＋ｘ［０］ のもとで解くわけだから 初期条件により ｘ［０］＝ａ と ｘ［１］＝ｂ・ｈ＋ａ が分かる （＊）においてｎ＝０として ｘ［２］＝２・ｘ［１］－（ｈ＾２＋１）・ｘ［０］＋ｈ＾２・ｕ［０］ は既に求まっているｘ［０］とｘ［１］から求まる （＊）においてｎ＝１として ｘ［３］＝２・ｘ［２］－（ｈ＾２＋１）・ｘ［１］＋ｈ＾２・ｕ［１］ は既に求まっているｘ［１］とｘ［２］から求まる （＊）においてｎ＝２として ｘ［４］＝２・ｘ［３］－（ｈ＾２＋１）・ｘ［２］＋ｈ＾２・ｕ［２］ は既に求まっているｘ［２］とｘ［３］から求まる （＊）においてｎ＝３として ｘ［５］＝２・ｘ［４］－（ｈ＾２＋１）・ｘ［３］＋ｈ＾２・ｕ［３］ は既に求まっているｘ［３］とｘ［４］から求まる （＊）においてｎ＝４として ｘ［６］＝２・ｘ［５］－（ｈ＾２＋１）・ｘ［４］＋ｈ＾２・ｕ［４］ は既に求まっているｘ［４］とｘ［５］から求まる ・・・・・・・・・・・・・・・・・ ｕ［ｎ］＝０（ｎ＝０，１，２，３，・・・）とすれば最初の問題が ｕ［ｎ］＝ｓｉｎ（ω・ｈ・ｎ） （ｎ＝０，１，２，３，・・・）とすれば後の問題が求まる
その他の回答 (4)
- nubou
- ベストアンサー率22% (116/506)
前回答に間違いが少なくとも一箇所あります 管理者に怒られるので指摘できません 補足で指摘してください
補足
ご回答ありがとうございます。 間違えというのはおそらく、 (x[n+1]-x[n]+x[n-1])/h^2+x[n]=0が (x[n+1]-2x[n]+x[n-1])/h^2+x[n]=0 ですよね？ 上の二つの場合で計算してみましたが、どちらも共振してしまいます。 なぜでしょうか？ 何度もすみませんがよろしくお願いします。
- nubou
- ベストアンサー率22% (116/506)
（ｘ［ｎ＋２］－ｘ［ｎ＋１］＋ｘ［ｎ］）／ｈ＾２＋ｘ［ｎ］＝ｕ［ｎ］ すなわち ｘ［ｎ＋２］＝ｘ［ｎ＋１］－（ｈ＾２＋１）・ｘ［ｎ］＋ｈ＾２・ｕ［ｎ］・・・（＊） を ｘ［０］＝ａ （ｘ［１］－ｘ［０］）／ｈ＝ｂすなわちｘ［１］＝ｂ・ｈ＋ｘ［０］ のもとで解くわけだから ｘ［０］＝ａ ｘ［１］＝ｂ・ｈ＋ａ （＊）においてｎ＝０として ｘ［２］＝ｘ［１］－（ｈ＾２＋１）・ｘ［０］＋ｈ＾２・ｕ［０］ ＝ｂ／ｈ＋ａ－（ｈ＾２＋１）・ａ＋ｈ＾２・ｕ［０］ （＊）においてｎ＝１として ｘ［３］＝ｘ［２］－（ｈ＾２＋１）・ｘ［１］＋ｈ＾２・ｕ［１］ （＊）においてｎ＝２として ｘ［４］＝ｘ［３］－（ｈ＾２＋１）・ｘ［２］＋ｈ＾２・ｕ［２］ （＊）においてｎ＝３として ｘ［５］＝ｘ［４］－（ｈ＾２＋１）・ｘ［３］＋ｈ＾２・ｕ［３］ （＊）においてｎ＝４として ｘ［６］＝ｘ［５］－（ｈ＾２＋１）・ｘ［４］＋ｈ＾２・ｕ［４］ ・・・・・・・・・・・・・・・・・ ｕ［ｎ］＝０（ｎ＝０，１，２，３，・・・）とすれば最初の問題が ｕ［ｎ］＝ｓｉｎ（ω・ｈ・ｎ） （ｎ＝０，１，２，３，・・・）とすれば後の問題が求まる なお式は保証できませんので自分で確かめてください 「自身あり」は考え方だけですので悪しからず
- nubou
- ベストアンサー率22% (116/506)
（ｘ［ｎ＋２］－ｘ［ｎ＋１］＋ｘ［ｎ］）／ｈ＾２＋ｘ［ｎ］＝ｕ［ｎ］ すなわち ｘ［ｎ＋２］＝ｘ［ｎ＋１］－（ｈ＾２＋１）・ｘ［ｎ］＋ｈ＾２・ｕ［ｎ］・・・（＊） を ｘ［０］＝ａ （ｘ［１］－ｘ［０］）／ｈ＝ｂすなわちｘ［１］＝ｂ・ｈ＋ｘ［０］ のもとで解くわけだから ｘ［０］＝ａ ｘ［１］＝ｂ／ｈ＋ａ （＊）においてｎ＝０として ｘ［２］＝ｘ［１］－（ｈ＾２＋１）・ｘ［０］＋ｈ＾２・ｕ［０］ ＝ｂ／ｈ＋ａ－（ｈ＾２＋１）・ａ＋ｈ＾２・ｕ［０］ 以下ｎ＝１，２，３，４，・・・とすればすべての０以上の整数ｎについて ｘ［ｎ］がもとまる ｕ［ｎ］＝０（ｎ＝０，１，２，３，・・・）とすれば最初の問題が ｕ［ｎ］＝ｓｉｎ（ω・ｈ・ｎ） （ｎ＝０，１，２，３，・・・）とすれば後の問題が求まる なお式は保証できませんので自分で確かめてください 「自身あり」は考え方だけですので悪しからず
- nubou
- ベストアンサー率22% (116/506)
ｘ（ｔ）→ｘ［ｎ］ ｘ’（ｔ）→ｘ１［ｎ］＝（ｘ［ｎ＋１］－ｘ［ｎ］）／ｈ ｘ”（ｔ）→ｘ２［ｎ］＝（ｘ１［ｎ＋１］－ｘ１［ｎ］）／ｈ ｕ（ｔ）→ｕ［ｎ］＝ｓｉｎ（ω・ｈ・ｎ） とおけばいいのでは？ すると ｘ”（ｔ）→ｘ２［ｎ］＝（ｘ［ｎ＋２］－ｘ［ｎ＋１］＋ｘ［ｎ］）／ｈ＾２ ｘ（０）＝ａ→ｘ［０］＝ａ ｘ’（０）＝ｂ→ｘ１［０］＝（ｘ［１］－ｘ［０］）／ｈ＝ｂ になります ｘ１［ｎ］とｘ２［ｎ］を消去してｘ［ｎ］だけの式にして解いてください ０＜ｈを小さくすればするほど近似はようなります
補足
ご回答ありがとうございます。 今頃遅いですが、自分で書いた差分方程式がメチャクチャでしたので訂正します。 (x[n+1]-2x[n]+x[n-1])/dt^2=-x[n] そして質問がございます。おそらくnubouさんが書いてくださった式を整理すると差分方程式は、 (x[n+2]-2x[n+1]+x[n])/h^2=-x[n] になりますよね？初期条件のx[0]=aは分かりますが、x1[0]=(x[1]-x[0])/h^2=bをどのように使用したらよいのか分かりません。また、x[1]はどのように求めるのでしょうか？連立方程式を解くのでしょうか？ お手数ですがよろしくお願いします。
関連するQ&A
- 波動方程式の差分法による境界条件
波動方程式(ρ*∂^2u/∂t^2=T*∂u^2/∂x^2)を差分法で、境界条件をx=0とx=Lで自由条件（T*∂u/∂x=0)とした場合を考えています。 上の波動方程式を差分化すると、 (u[n+1][j]-2*u[n][j]+u[n-1][j])/(dt^2)=(u[n][j+1]-2*u[n][j]+u[n][j-1])/(dx^2) の形になると思います。（T/ρ=1.0とし、nは時間、jは距離の格子点として考えています） 初期条件は適当な形状を与えます。 両端を固定条件（u[n][0]=0,u[n][xp]=0,）とした場合はうまく解を得ることが出来ました。（xpはx=Lでの格子点） 問題は自由条件（T*∂u/∂x=0)すなわち、 T*(u[n][1]-u[n][0])/dx=0、T*(u[n][xp+1]-u[n][xp])/dx=0 となる場合、これをどのように使用したらよいのでしょうか？ または根本的に考え方が間違っているのでしょうか？ 本当に困ってます。よろしくお願いいたします。 内容が不十分の場合は補足要求お願いします。
- ベストアンサー
- その他（学問・教育）
- 差分法
過去にも同様な質問があったようですが、それでも理解できないのでご教授お願いします。 拡散方程式などで 境界条件が：∂u/∂x=0 となる場合、 前進差分で表すと：(u[n][j+1]-u[n][j])/δx=0 後退差分で表すと：(u[n][j]-u[n][j-1])/δx=0 中心差分で表すと：(u[n][j+1]-u[n][j-1])/2*δx=0 となります。（n,jは時間、位置の格子点） それぞれ整頓すると、j=0（左の境界点）では 前進差分：u[n][1]=u[n][0] 後退差分：u[n][0]=u[n][-1] 中心差分：u[n][1]=u[n][-1] となるのですがぁ、 本などではよく中心差分を使っているのですがこの理由がよく分かりません。 上の3つの式から：u[n][-1]=u[n][0]=u[n][1]としてしまってもよいのでしょうか？ 以上2点と何かアドバイスお願いします。
- ベストアンサー
- その他（学問・教育）
- 偏微分方程式の問題です。
・問題 u=u(t,x) u_t=ku_｛xx｝ k＞0 (0＜x＜1,t＞0) 初期条件 u(0,x) sin(πx)+1/2sin(3πx) 境界条件 u(t,0)=u(t,1)=0 ・答え 境界条件より正弦級数展開をする． （☆）u(t,x)=Σ_{n=1}^∞b_n(t)sin(nπx) （★）b_n(t)=2∫_0^1u(t,x)sin(nπx)dx すると， u_t=Σ_{n=1}^∞{db_n(t)/dt}sin(nπx) ku_{xx}=Σ_{n=1}^∞b_n(t)d^2{sin(nπx)}/dx^2 =Σ_{n=1}^∞b_n(t)k(-nπ)^2sin(nπx) u_t=ku_{xx}より db_n(t)/dt=b_n(t)k(-nπ)^2=-n^2π^2kb_n(t) ∴b_n(t)=b_n(0)e^{-n^2π^2kt} ☆に代入して （☆☆）u(t,x)=Σ_{n=1}^∞b_n(0)e^{-n^2π^2kt}sin(nπx) ∴u(0,x)=Σ_{n=1}^∞b_n(0)sin(nπx) これと初期条件u(0,x)=sin(πx)+(1/2)sin(3πx)を係数比較して b_1(0)=1,b_3(0)=1/2,b_n(0)=0(n≠1,3) ☆☆に代入して ∴u(t,x)=e^{-π^2kt}sin(πx)+(1/2)e^{-9π^2kt}sin(3πx)（答） ーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーーー 以前、上記の問題をここで質問し、解答を頂いたのですが、この部分の途中式がどうしてもわかりません。 ku_{xx}=Σ_{n=1}^∞b_n(t)d^2{sin(nπx)}/dx^2 =Σ_{n=1}^∞b_n(t)k(-nπ)^2sin(nπx) お手数ですが宜しくお願い致します！
- ベストアンサー
- 数学・算数
- 偏微分方程式、差分法
Fitzhugh-Nagumo方程式 dU/dt = d/dx(dU/dx) + (a-U)(U-1)U-V (d/dx(dU/dx)はUのxに対する二回微分) dV/dt =e(bU-gV) a=0.1, b=0.5, g=1.0, e=0.01 初期条件 (U,V)=(1,0) if x=0 (U,V)=(0,0) if x>0 境界条件 dU/dx=0 at x=0 dV/dx=0 at x=0 を差分して陽解法で解くと添付の答えになるらしいのですが、自分で解いたところ添付の結果のようになりました。答えと一致しないため、プログラム上で何を間違えているのかご指摘頂けると助かります。 #include <stdio.h> #include <stdlib.h> //exsit()で必要 #include <math.h> int main(){ double a,b,g,e; double dt,dx; int x,t; double **u,**v; int i,j,k; a=0.1;b=0.5;g=1.0;e=0.01; x=100;t=5000.; dt=0.001;dx=0.1; FILE *fp; if((fp = fopen("ResultNagumo.txt","w"))==NULL){ printf("Can't open file\n"); exit(2); } u = (double**)malloc(sizeof(double*)*t); v = (double**)malloc(sizeof(double*)*t); for(i=0;i<t;i++){ u[i]=(double*)malloc(sizeof(double)*x); v[i]=(double*)malloc(sizeof(double)*x); } //初期値 u[0][0]=1.0; v[0][0]=0.0; for(i=1;i<x;i++){ u[0][i]=0.0; v[0][i]=0.0; } //差分計算 for(i=0;i<t-1;i++){ for(j=1;j<x-1;j++){ u[i+1][j] = u[i][j] + dt*((u[i][j+1]-2*u[i][j]+u[i][j-1])/pow(dx,2)+(a-u[i][j])*(u[i][j]-1)*u[i][j]-v[i][j]); v[i+1][j] = v[i][j] + dt*(e*(b*u[i][j]-g*v[i][j])); } //境界条件 u[i+1][0]=u[i+1][1]; v[i+1][0]=v[i+1][1]; u[i+1][x-1]=u[i+1][x-2]; v[i+1][x-1]=v[i+1][x-2]; } //結果の出力 for(i=0;i<t;i++){ if((i%100)==0){ fprintf(fp,"%d\n",i); for(j=0;j<x;j++){ fprintf(fp,"%2.4e,",u[i][j]); } } } for(i=0;i<t;i++){ free(u[i]); free(v[i]); } free(u);free(v); fclose(fp); getchar(); return 0; }
- ベストアンサー
- C・C++・C#
- 積分すると1/2とLが消える？
本の計算では、ある式を積分すると1/2とLが消えてしまっていて、自分の計算と合いません。右辺が0なので両辺に2を掛けてLで割った可能性もあるのですが、1/2が消えるのはまだ許せても、Lが消えるのは納得いきません。（消していいものなのか、もしくは別の理由で消えたのか、判断願います。） 原文を引用します： ※関係式： ∫[0,L] cos (2nπx/L) dx = ∫[0,L] sin (2nπx/L) dx = 0 (式1.17a) (nは正あるいは負の整数) [前略] …さらに、 (1/2) * { da[0](t) }/dt + Σ[n=1,∞] { da[n](t) }/dt * cos(2nπx/L) + Σ[n=1,∞] { db[n](t) }/dt * sin(2nπx/L) = -D * Σ[n=1,∞] {(2nπ/L)^2} * a[n](t) * cos(2nπx/L) - D * Σ[n=1,∞] {(2nπ/L)^2} * b[n](t) * sin(2nπx/L) （式7.6） …の両辺をそのままx=0からx=Lまで積分すると、先の関係式(式1.17a)より、 { da[0](t) }/dt = 0 (式7.7c) が得られる。 ・・・以上、引用終わり。 私の計算だと、（式7.6）の両辺をそのままx=0からx=Lまで積分するので： (1/2) * { da[0](t) }/dt *∫[0,L] 1 dx + Σ[n=1,∞] { da[n](t) }/dt * ∫[0,L] cos(2nπx/L) dx + Σ[n=1,∞] { db[n](t) }/dt * ∫[0,L] sin(2nπx/L) dx = -D * Σ[n=1,∞] {(2nπ/L)^2} * a[n](t) * ∫[0,L] cos(2nπx/L) dx - D * Σ[n=1,∞] {(2nπ/L)^2} * b[n](t) * ∫[0,L] sin(2nπx/L) dx になり、 ∫[0,L] cos (2nπx/L) dx = ∫[0,L] sin (2nπx/L) dx = 0 (式1.17a) により、 ∫[0,L] cos(2nπx/L) dx ∫[0,L] sin(2nπx/L) dx はすべて0になります。 残りを計算すると (1/2) * { da[0](t) }/dt *∫[0,L] 1 dx = 0 (1/2) * { da[0](t) }/dt * [x][0,L] = 0 (1/2) * { da[0](t) }/dt * [L-0] = 0 (1/2) * { da[0](t) }/dt * L = 0 になります。本の答えは { da[0](t) }/dt = 0 なので合いません。 [概略]（原文が凄まじく長いので、自分の言葉で書きました…なんとなく分かっていただけたらと思います…足りなかったら補足します） u(x,t) が時刻tにおける座標xでの温度uを表します。 このとき、1次元の熱伝導方程式 { δu(x,t) }/δt = D * { δ^(2) * u(x,t) }/{ δx^(2) } （式7.1） を解くためにフーリエ級数をどう使うのか、というのが今回のテーマです。 最初の時刻(t=0)での温度分布（つまり、初期条件）は u(x,0) = f(x) （式7.2） とします。 例として、長さLの「リング状の」棒での熱伝導を考えます。リング状なので、棒に沿って「ある点」をx=0とすると、 x=x[0] と x=x[0]+L は同じ点に対応するため、「温度uが座標xについて周期Lの周期関数である」という周期的境界条件 u(x,t) = u(x+L,t)（式7.3） が任意の時刻t (>=0)で成り立っている必要があります。 このときには初期温度分布f(x)にも条件が付き、 f(x) = f(x+L)（式7.4） が成り立っていないといけません。すなわち、f(x)も周期Lの周期関数です。 u(x,t)をフーリエ級数展開すると、 u(x,t) = (1/2) * a[0](t) + Σ[n=1,∞] a[n](t) cos (2nπx/L) + Σ[n=1,∞] b[n](t) sin (2nπx/L)（式7.5） になり、この（式7.5）を（式7.1）に代入して項別の偏微分をすると、私の質問に出てくる (1/2) * { da[0](t) }/dt + Σ[n=1,∞] { da[n](t) }/dt * cos(2nπx/L) + Σ[n=1,∞] { db[n](t) }/dt * sin(2nπx/L) = -D * Σ[n=1,∞] {(2nπ/L)^2} * a[n](t) * cos(2nπx/L) - D * Σ[n=1,∞] {(2nπ/L)^2} * b[n](t) * sin(2nπx/L) （式7.6） になります。
- ベストアンサー
- 数学・算数
- 線形差分方程式の一般解について
2階の線形差分方程式 x(n+2)-a*x(n+1)-b*x(n)=0 は特性方程式の２解をαとβとおくと、 α≠βのとき x(n)=A1*α^n+A2*β^n （１） となり重解のとき、 x(n)=A1*α^n+A2*n*α^n （２） となりますが、両式の違いであるnは どのように導かれたのでしょうか。 またA1とA2の求め方も教えてもらえると助かります。 お願いします。
- ベストアンサー
- 数学・算数
- 減衰振動の微分方程式の解
先生から配られたプリントには減衰振動の微分方程式が「m(dx/dt)^2+2γ・dx/dt+ω^2x=0」の時、解が「x=A・EXP（-γt）cos（ω´t＋φ）」って書かれてます。 摩擦：гならг/m=2γ、バネ定数：kならk/m=ω^2、A=√C1^2+C2^2、φ=C1/C2、ω´^2=ω^2-γ^2です。 解の式で、cosじゃなくてsinではないのですか？単振動・強制振動の場合も同様にcosでした。 誰かよろしくお願いいたします。
- ベストアンサー
- 物理学
- 拡散方程式の数値計算(有限差分法)による解き方
下記の拡散方程式を数値計算(有限差分法)で解き方について教えてください。 拡散方程式 D∂^2C(x,y,t)/∂x^2+u∂C(x,y,t)/∂x=∂C(x,y,t)/∂t 境界条件 D∂C(x,y,t)/∂x=(K/β-u)×C(x,y,t)-KC'(y,t) at x=0 初期条件 C(x,y,0)=Co at t=0 質量保存則 ∂C'(y,t)/∂y=4/vd(K/β×C(0,y,t)-KC'(y,t) ) ---------------------------------------------------------- また、可能であれば、有限差分法以外にも数値解法(有限要素法、有限体積法など)があると思いますが、拡散方程式を解く場合、どの解法が一般的に適しているのか教えてください！ よろしくお願い申し上げます。
- 締切済み
- 数学・算数
- 物理学の単振動の問題を教えてください
この物理の単振動の問題を教えてください。 「以下の写真のように、人と、水平方向に単振動をする球体があるとします。 人がこれを押す外力の周期が球体の周期T(単振動ω)と同じとき、球体の運動を考え以下の問いに答えなさい。ただし、ω=1と考えて答えなさい。 (1)単位質量あたりの外力をf(t)として鉛直真下を原点とした変位xに対する運動方程式(微分方程式)を書きなさい。 ➁外力f(t)が下の写真のグラフの時、0≦t≦3Tの解を求めなさい。ただし、t=0の時x=-1、dx/dt=0とする。また、t=T/2、T、3T/2T、5T/2、3Tでは、x、dx/dtが連続とする。 ➂0≦ｔ≦3Tの解のグラフを書きなさい。」 一部でも構いません。分かる方、教えてください
- ベストアンサー
- 物理学
補足
ご回答ありがとうございます。遅くなりましてすみません。 更正された式で試してみませたが、計算結果は以前の式のものと全く同じになりました。難しいですねぇ。。まだ何か方法がありましたらお願いいたします。