常微分方程式

常微分方程式 傾きがわかった関数の軌跡を求める．

微分方程式の種類 • 常微分方程式→独立変数が一個単振動 m:mass, c:dash pot, k:spring, t:time, x:axis • 偏微分方程式→独立変数が２個以上　　　　　　　　　　　　　　　　　　　一次元熱伝達 T:temperature, t:time, x:axis 二階の常微分方程式 d2x dx m + c + kx = 0 dt2 dt 二階の偏微分方程式 ∂T ∂2T cr = l ∂t∂x2

微分方程式を解くという意味 dy dx わかっている各点の傾きが既知 y この線を求める初期値初期値 x

微分方程式を解くという意味 dy = f(x,y) dx 与えられる式与えられる初期条件 y(x0) = y0 求める値 y(x)

オイラー法 y2 = y1 + y’(x1,y1)h y y1 = y0 + y’(x0,y0)h h x y’(x1,y1)h 初期値 (x0,y0) y’(x0,y0)h x1 x3 x2 傾きy’に，区間幅hをかけて，たす．

オイラー法のアルゴリズム Start h = (xn-x0)/n yi+1 = yi + y’(xi)*h No xi = xn? Yes End

テキストの例題 初期値(0,0) = 2 x dy dx 傾きがXに比例 y この線を求める初期値 x 10 0

Eulerの公式による積分 #include <stdio.h> #include <math.h> #define EPS .00000001 double func_f(double); int main( ){ double x=0.0, y=0.0; double h=0.01, dx=1.0, xmax=10.0; double ddx=0.0; printf(" X\t Y\n");

do{ if(x>=ddx-EPS){ ddx+=dx; printf("%7.4lf %7.4lf\n", x, y); } y+=h*func_f(x); x+=h; }while(x<=xmax); return 0; } double func_f(double x) { return 2.0*x; } １おきに出力 yi+1 = yi + y’(xi)*h xi = xn? y = x2 y’ = 2x

計算結果

do{ if(x>=ddx-EPS){ ddx+=dx; printf("%7.4lf %7.4lf\n", x, y); } y+=h*func_f(x); x+=h; }while(x<=xmax); return 0; } double func_f(double x) { return 2.0*x; } １おきに出力不要な部分を取って計算する

計算結果 X 0.0000 0.0100 0.0200 10.0000 Y 0.0000 0.0000 0.0006 99.9000 dy dx = 2 x 積分 y= x2 y = 102 = 100 同じはず×

Ｅｘｃｅｌで確認 dy = 2 x dx 区間の幅前の値＋傾き×幅

Ｅｘｃｅｌで確認 dy = 2 x dx 区間の幅前の値＋傾き×幅プロット

解析解との比較

解析解との比較 y 100 y = x2 誤差 75 50 25 x 10 0 5

オイラー法 次の値は，前の値に、傾きと区間幅をかけたものを足す．精度は分割数を増やすと高くなるす．

ｘ，ｙ変数をベクトルに変更 double x[1001], y[1001]; double h, xmax=10.0; int i, n=1000; x[0]=0.0; y[0]=0.0; h = (xmax-x[0])/((double) n); printf(" X\t Y\n"); for(i=0;i<n;i++){ y[i+1] = y[i]+h*func_f(x[i]); x[i+1] = x[i]+h; } for(i=0;i<=n;i+=100){ printf("%7.4lf %7.4lf\n", x[i], y[i]); }

力学への応用（単振動） エクセルで解く • 単振動質量：m = 5 kg ばね定数：k = ○○○ N/m • 初期条件：x0 =0.1 m, v0 = 0.0 m/s k d2x F = m = - kx dt2 -kx m m m m F = ma d2x kx = - dt2 m x 周期T?

Ｅｘｃｅｌで計算 前の値＋傾き×幅周期T ?

Ｅｘｃｅｌで計算 前の値＋傾き×幅周期T ? プロット

波形が発散する 刻み幅を縮小する波形が発散する理由は何か？誤差 h = 0.001 で再計算！

減衰比A1/A0 ? 周期T ? 減衰比A1/A0 A0 A1 周期T

アルゴリズム 来週の課題Ｐ７３ Start ai = -(cvi + kxi)/m vi+1 = vi + ai*h xi+1 = xi + vi*h Yes t < tmax? No End

力学への応用 d2x F = - μ m g= m dt2 d2x = - μ g dt2 x F 自動車の停止するまでの距離質量：m = 2000 kg 摩擦係数：μ = 0.8 (g = 9.8 m/s2) 初期条件：v0 = 25 m/s, x0 = 0 m

アルゴリズム Start vi+1 = vi + ai*h xi+1 = xi + vi*h Yes vi+1 > 0 ? No End

Ｅｘｃｅｌで計算 前の値＋傾き×幅

Ｅｘｃｅｌで計算 前の値＋傾き×幅プロット

Cプログラム #include <stdio.h> int main( ){ double x=0.0, v=25.0, h=0.01, t=0.0; double g=9.8, mu=0.8; printf(" T\t X\t V\n"); do{ printf(“ %7.4lf\t %7.4lf %7.4lf\n”, t, x, v; t=t+h; v=v-mu*g*h; x=x+v*h; }while(v>0.0); return 0; }

常微分方程式

常微分方程式

Presentation Transcript