常微分方程式と偏微分方程式 １．常微分方程式 独立変数が一個のもの 振動の運動方程式 ２．偏微分方程式 独立変数が二個以上のもの 　　独立変数が一個のもの 二階の常微分方程式 振動の運動方程式 ２．偏微分方程式 　　独立変数が二個以上のもの 二階の偏微分方程式 一次元熱伝導

一階の常微分方程式の解析 初期値から解を求める。 与えられる式 与えられる初期条件 求める値

オイラー法 h = xn/n y yn x1 x y0 xn x0 y1 = y0 + y’(x0,y0)h 初期値が既知なので、逐次的に値を求める。 h = xn/n y yn x1 x y0 xn x0 y1 = y0 + y’(x0,y0)h

オイラー法の流れ x0=0, y0=1 y(1)? START 例題：dy/dx=y+x2-2x h=(xn-x0)/n yi+1=yi+yi’h i=n END No Yes 例題：dy/dx=y+x2-2x x0=0, y0=1 y(1)? y x0 y0 x1 h x xn y1 = y0 + y’(x0,y0)h

エクセルで作ってみよう（１） 時間： 位置： 空気密度： 投影面積： 抗力係数： 水滴の落下

エクセルで作ってみよう（１） 水滴の落下 水滴半径： 水密度： 球とすると とすると

=-$B$6*$B$5*$B$4/$B$3*E2^2+$B$8 A B C D E F 1 半径 0.002 時刻 位置 速度 加速度 2 水密度 μ 1000 3 質量 m 4 投影面積S 5 抗力係数CD 0.4 6 空気密度ρ 1.4 7 時間刻みΔt 0.01 8 重力加速度g 9.8 =-$B$6*$B$5*$B$4/$B$3*E2^2+$B$8 コピー =C2+$B$7 =4/3*PI()*B1^3*B2 =PI()*B1^2

=D2+0.5*(E2+E3)*$B$7 =E2+F2*$B$7 Dt dx dt d2x dt2 A B C D E F 1 半径 0.002 時刻 位置 速度 加速度 2 水 密度 1000 9.8 3 質量 3.35E-5 0.01 4 投影面積 1.26E-5 0.02 5 抗力係数 0.4 0.03 6 空気密度 1.4 0.04 7 時間刻み 0.05 8 重力加速度 0.06 微分値から次の値を求める =D2+0.5*(E2+E3)*$B$7 =E2+F2*$B$7 Dt dx dt d2x dt2

Cプログラムへ A B C D E F 1 半径 0.002 時刻 位置 速度 加速度 2 水 密度 1000 9.8 3 質量 9.8 3 質量 3.35E-5 0.01 4 投影面積 1.26E-5 0.02 5 抗力係数 0.4 0.03 6 空気密度 1.4 0.04 7 時間刻み 0.05 8 重力加速度 0.06 すべて配列にする すべて変数にする

int i,n=200; double r=0.002, dw=1000.0, da=1.4; double dt=0.01, m, s, cd=0.4, g=9.8 double t[200], x[100], v[200], a[200]; m=4.0/3.0*3.1415*r*r*r; s=3.1415*r*r; t[0]=0.0; x[0]=0.0; v[0]=0.0; a[0]=g; for(i=1;i<n;i++){ t[i]=t[i-1]+dt; v[i]=v[i-1]+a[i-1]*dt; x[i]=x[i-1]+v[i-1]*dt; a[i]=-cd*da*s/m*v[i]*v[i]+g;

振動のシミュレーション １．常微分方程式 k =100 N/m c = 500N/(m/s) t = 0 (s)のとき m = 2 kg x = 0.2 m v = 0 m/s m = 2 kg m

振動のシミュレーション １．常微分方程式 k =100 N/m c = 500N/(m/s) t = 0 (s)のとき m = 2 kg x = 0.2 m v = 0 m/s m = 2 kg

=-1/$B$3*($B$2*E2+$B$1*D2) A B C D E F 1 ばね定数 k 100 時刻 位置 速度 加速度 2 粘性係数 c 500 3 質量 m 4 時間刻み 0.01 5 初期位置 0.0 6 初期速度 7 =-1/$B$3*($B$2*E2+$B$1*D2) =A5 =A6 コピー =C2+$B$4

FILE *fp = fopen("data.csv","w"); int i,n=20000; double k=100.0, c=500.0, dt=0.01, m=2.0; double t[20000], x[20000], v[20000], a[20000]; t[0]=0.0; x[0]=0.2; v[0]=0.0; a[0]= -1.0/m*(c*v[0]+k*x[0]); for(i=1;i<n;i++){ t[i]=t[i-1]+dt; v[i]=v[i-1]+a[i-1]*dt; x[i]=x[i-1]+v[i-1]*dt; a[i]=-1.0/m*(c*v[i]+k*x[i]); fprintf(fp,"%d,%f,%f\n",i,t[i],x[i]); } fclose(fp);

熱伝導の基礎 温度T x x l:熱伝導率 A:断面積 伝わる熱量は温度勾配に比例

偏微分方程式 例：熱伝導 x x l:熱伝導率 c:比熱 r:密度 A:断面積 流入する熱量は温度上昇を引き起こす。

差分法による解法 Dt:時間における熱の移動 T(i,j - 1) T(i,j) T(i,j + 1) x h h

差分法による解法 Dt:時間における熱の移動 T(i,j - 1) T(i,j) T(i,j + 1) x h h h h 前と後ろの温度を足して、現在の温度の２倍を引く。

t = 0, x<1m;T=0Co, x=1;T=100Co 偏微分方程式 T = 100 Co x 1 t = 0, x<1m;T=0Co, x=1;T=100Co

=D3+$B$1/($B$2*$B$3)*(E3-2*D3+C3)*$B$4/$B$7^2 A B C D E F N 1 熱伝導率 W/mK 390 座標 2 比熱 J/kgK 380 時刻 0.1 0.2 1.0 3 密度 Kg/m3 8900 4 時間刻み 0.01 5 初期温度 ℃ 0.0 6 端点温度 100 7 X座標 刻み =D2+$B$7 =$B$5 =$B$6 =C3+$B$4 =D3+$B$1/($B$2*$B$3)*(E3-2*D3+C3)*$B$4/$B$7^2

t[i][j], 時間 位置 時間の数 要素の数 int i, j, m=100, n=10; double x[100], t[1000][100], h, dt=5.0; double k=390.0, c=380.0, r=8900.0, t1=100.0; x[0]=0.0, x[n]=1.0; for(i=0,i<m+1;i++){ t[i][0]=0.0; t[i][n]=t1; } h = (x[n] - x[0])/(double)n; for(j=0,j<n;n++){ x[j+1] = h * (double) (j + 1) t[0][j]=0.0; for(i=0;i<m;i++){ for(j=1;j<n;j++){ t[i+1][j] = t[i][j] + k/(c*r)*(t[i][j-1]-2.0*t[i][j]+t[i][j+1])/(h*h)*dt; return 0; 温度の変数 t[i][j], 時間の変数 要素の変数