傾きがわかった関数の軌跡を求める． 変数は二つ以上 偏微分方程式 傾きがわかった関数の軌跡を求める． 変数は二つ以上

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

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

変数が２個の場合 ∂u ∂x ，　　　　わかっている ∂y u 初期値 y 傾きが分っている この面を求める x

変数が２個の場合 ∂u ∂x ，　　　　わかっている ∂y u 初期値 y Dx 次の点を 求める方法 x

変数が２個の場合 ∂u u ui - ui-1 = Δx ∂x x = xi x = xi-1 y yj-2 yj-1 yj yj+1 次の点 x = xi x = xi-1 前の点 y yj-2 yj-1 yj yj+1 yj+2 yj+3

二つの変数の例 （偏微分方程式の前に） h1, h2, h3, h4, h5の変化？ 初めは，すべて h1～h5 = 1.0 時間に対する 高さの変化？ h1 h2 h3 h4 h5 高さh 単位時間 当たりの流量 F = h1 – h2 F = h1 – h2

二つの変数の例 （偏微分方程式の前に） 前の点 h1 h2 h3 h4 h5 高さh 次の点

二つの変数の例 （偏微分方程式の前に） h1, h2, h3, h4, h5の変化？ 初めは，すべて h1～h5 = 1.0 時間に対する 高さの変化？ h1 h2 h3 h4 h5 高さh Dt 時間の増加量 Dt (h2-h3)- Dt (h3-h4) 断面積 A = 0.5 単位時間 当たりの流量 F = h1 – h2 F = h1 – h2

二つの変数の例 （偏微分方程式の前に） h1 h2 h3 h4 h5 高さh 単位時間 当たりの流量 F = h4 – h5

二つの変数の例 （偏微分方程式の前に） h1 h2 h3 h4 h5 高さh 単位時間 当たりの流量 F = h4 – h5

二つの変数の例 （偏微分方程式の前に） h1 h2 h3 h4 h5 高さh 単位時間 当たりの流量 F = h4 – h5

二つの変数の例 （偏微分方程式の前に） h1 h2 h3 h4 h5 高さh 単位時間 当たりの流量 F = h4 – h5

二つの変数の例 （偏微分方程式の前に） h1 h2 h3 h4 h5 高さh 単位時間 当たりの流量 F = h4 – h5

二つの変数の例 （偏微分方程式の前に） h1 h2 h3 h4 h5 高さh 単位時間 当たりの流量 F = h4 – h5

二つの変数の例 （偏微分方程式の前に） h1, h2, h3, h4, h5の変化？ 初めは，すべて h1～h5 = 1.0 時間に対する 高さの変化？ h1 h2 h3 h4 h5 高さh Dt 時間の増加量 Dt (h2-h3)- Dt (h3-h4) 断面積 A = 0.5 単位時間 当たりの流量 F = h1 – h2 F = h1 – h2

Ｅｘｃｅｌで計算 初期値 時間刻みDt A B C D E F 1 T h1 h2 h3 h4 h5 2 3 0.01 = B2 - 3 0.01 = B2 - = D2 – + 4 0.02 5 0.03 6 0.04 7 0.05 初期値 流出量/断面積 流入量/A 流出量/A 時間刻みDt

Ｅｘｃｅｌで確認 F = h1 – h2 Dt Dt A A B C 1 T h1 h2 2 3 0.01 3 0.01 =B2 - 0.01 *SQRT(B2 - C2) / 0.5 4 0.02 5 0.03 6 0.04 7 0.05 8 0.06 Dt Dt A 前の値 ー 流出 前の値 ー 流出 ＋ 流入

Ｅｘｃｅｌで確認 F = h1 – h2 Dt Dt A A B C 1 T h1 h2 2 3 0.01 =B2 - 0.01 3 0.01 =B2 - 0.01 *SQRT(B2 - C2) / 0.5 =C2 - 0.01*SQRT(C2 - D2) / 0.5 +0.01*SQRT(B2 – C2) / 0.5 4 0.02 5 0.03 6 0.04 Dt Dt A 前の値 ー流出 前の値 ー流出＋流入

計算結果 1秒後 2秒後

二つの変数の例 （偏微分方程式の前に） 一秒間に １ｍ３ h1, h2, h3, h4, h5の変化？ 初めは，すべて 十分時間が たったときの h1～h5　？ h1 h2 h3 h4 h5 高さh 断面積 A = 0.5 単位時間 当たりの流量 F = h1 – h2

Ｅｘｃｅｌで計算 初期値 A B C D E F 1 T h1 h2 h3 h4 h5 2 3 0.01 =B2 - 0.01 3 0.01 =B2 - 0.01 *SQRT(B2 - C2) / 0.5 + 4 0.02 5 0.03 6 0.04 7 0.05 初期値 時間刻み当たりの流入高さ 前回と同じ

Ｅｘｃｅｌで計算 初期値 無限 A B C D E F 1 T h1 h2 h3 h4 h5 2 3 0.01 =B2 - 0.01 3 0.01 =B2 - 0.01 *SQRT(B2 - C2) / 0.5 + 4 0.02 5 0.03 6 0.04 7 0.05 無限 H1? H2? H3? H4? H5? 初期値 時間刻み当たりの流入高さ 今日の課題 前回と同じ

Ｅｘｃｅｌで計算 初期値 A B C D 1 T h1 h2 h3 2 3 0.01 =B2 - 0.01 3 0.01 =B2 - 0.01 *SQRT(B2 - C2) / 0.5 +0.01*1/0.5 =C2 - 0.01*SQRT(C2 - D2) / 0.5 +0.01*SQRT(B2 – C2) / 0.5 4 0.02 5 0.03 6 0.04 7 0.05 時間刻み 当たりの 流出量/断面積 前回と同じ

Ｅｘｃｅｌで計算 答え A B C D E F 1 T h1 h2 h3 h4 h5 2 3 0.01 0.02 - 10 3.432258 3 0.01 0.02 - 10 3.432258 2.571056 1.826114 1.169101 0.57019 100 4.999046 3.999123 2.999272 1.999479 0.999729 無限 5 4 答え

熱伝導に応用 T1 T2 T3 T4 T5 高さh 温度T h Q = l A(T4 – T5)/h x l:熱伝導率 c:比熱 r:密度 単位時間 当たりの熱の流れ Q = l A(T4 – T5)/h

一区間の温度上昇 T1 T2 T3 T4 T5 T4’=T4+ T ÷ρcAh h h 熱容量:=ρcAh +l A(T3 – T4)/h x 熱容量:=ρcAh T1 T2 T3 T4 T5 T4’=T4+ +l A(T3 – T4)/h T -l A(T4 – T5)/h ÷ρcAh h h l:熱伝導率 c:比熱 r:密度 A:断面積 Q = l A(T4 – T5)/h Q = l A(T3 – T4)/h

1 T4’=T4+ w1 w2 w3 w4 w5 u1 u2 u3 u4 u5 ÷ρcAh T w[i]=u[i]+ h h +l A(T3 – T4)/h -l A(T4 – T5)/h u1 u2 u3 u4 u5 ÷ρcAh T w[i]=u[i]+ +(u[i-1]–u[i])/h2 h h -(u[i]-u[I+1])/h2 Q = l A(u3 – u4)/h Q = l A(u4 – u5)/h

#include <stdio.h> #define N 20 int main( ) { double u[N+1], w[N+1]; double k=0.001; double h, r, s; int i, j; h = 1.0/(double)N; r = k/(h*h); s = 1.0 - 2.0*r; 今の温度 次の温度 時間刻み w[j] w, u t = ti u[j] u[j+1] ｘの区間 t = ti-1 u[j-1] x xj-2 xj-1 xj xj+1 xj+2 xj+3

l ui+1 – 2ui + ui-1 cr wi = ui + Dt h2 for(i=0;i< N;i++) u[i]=1.0; for(i=0;i<=N;i++) w[i]=0.0; u[0] = 0.0; u[N] = 0.0; for(j=1;j<=200;j++){ if((j%10)==0){ printf("%5.3lf ", (double)j*k); for(i=0;i<=N;i+=2) printf("%5.3lf ", u[i]); printf("\n"); } for(i=1;i<N;i++) w[i] = r*(u[i+1]+u[i-1]) + s*u[i]; for(i=1;i<=N;i++) u[i] = w[i]; return 0; 初期条件 温度＝１（両端以外） 温度＝０（両端） 剰余演算子 １０で割ったあまり 出力 r = k/(h*h); s = 1.0 - 2.0*r; 1 wi = ui + Dt ui+1 – 2ui + ui-1 h2 l cr

for(i=1;i<N;i++) w[i] = r*(u[i+1]+u[i-1]) + s*u[i]; for(i=1;i<=N;i++) u[i] = w[i]; h = 1.0/N r = k/(h*h) s = 1.0 - 2.0*r = 1 - 2*k/h2 ui+1 – 2ui + ui-1 h2 wi = ui + Dt

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

熱伝導 ∂T ∂x ∂T ∂x Q’ = - l A + (- l A) Dx ∂ Q = - l A x + Δx Q - Q’ = l ADx = cr ADx ∂2T ∂x2 ∂T ∂t l:熱伝導率 c:比熱 r:密度 A:断面積 体積 = ∂2T ∂x2 ∂T ∂t l cr xに関する二階微分

の式 T t = ti t = ti-1 h x ∂2T ∂x2 Tj+1 - Tj Tj - Tj-1 - h h h xj-2 xj-1 x+Δx T t = ti t = ti-1 前の点 ここの２階微分 Tj+1 - Tj Tj - Tj-1 - h h h h x xj-2 xj-1 xj xj+1 xj+2 xj+3 iは時間変化 jは座標変化 ∂2T ∂x2 Tj+1 - 2Tj + Tj-1 = h2

Tiの式 ∂T ∂t Ti - Ti-1 = Dt = ∂2T ∂x2 ∂T ∂t l cr ∂2T ∂x2 Tj+1 - 2Tj + Tj-1 = h2 Ti = Ti-1 + Dt Tj+1 - 2Tj + Tj-1 h2 l cr

Tiの式 Ti T t = ti t = ti-1 Ti-1 x Ti = Ti-1 + Dt Tj+1 - 2Tj + Tj-1 h2 l cr Ti T t = ti t = ti-1 Ti-1 前の点 x xj-2 xj-1 xj xj+1 xj+2 xj+3

テキスト問題 t = 0 1 温度u w 1 l:熱伝導率 A:断面積 l cr = 1 両端の温度は常にゼロ

#include <stdio.h> #define N 20 int main( ) { double u[N+1], w[N+1]; double k=0.001; double h, r, s; int i, j; h = 1.0/(double)N; r = k/(h*h); s = 1.0 - 2.0*r; 今の温度 次の温度 時間刻み w[j] w, u t = ti u[j] u[j+1] ｘの区間 t = ti-1 u[j-1] x xj-2 xj-1 xj xj+1 xj+2 xj+3

l Tj+1 - 2Tj + Tj-1 cr Ti = Ti-1 + Dt h2 for(i=0;i< N;i++) u[i]=1.0; for(i=0;i<=N;i++) w[i]=0.0; u[0] = 0.0; u[N] = 0.0; for(j=1;j<=200;j++){ if((j%10)==0){ printf("%5.3lf ", (double)j*k); for(i=0;i<=N;i+=2) printf("%5.3lf ", u[i]); printf("\n"); } for(i=1;i<N;i++) w[i] = r*(u[i+1]+u[i-1]) + s*u[i]; for(i=1;i<=N;i++) u[i] = w[i]; return 0; 初期条件 温度＝１（両端以外） 温度＝０（両端） 剰余演算子 １０で割ったあまり 出力 r = k/(h*h); s = 1.0 - 2.0*r; 1 Ti = Ti-1 + Dt Tj+1 - 2Tj + Tj-1 h2 l cr

for(i=1;i<N;i++) w[i] = r*(u[i+1]+u[i-1]) + s*u[i]; for(i=1;i<=N;i++) u[i] = w[i]; h = 1.0/N r = k/(h*h) s = 1.0 - 2.0*r = 1 - 2*k/h2 Tj+1 - 2Tj + Tj-1 h2 Ti = Ti-1 + Dt

微分方程式を考えない方法 t = 0 1 温度u w 1 l:熱伝導率 A:断面積 l cr = 1 両端の温度は常にゼロ

微分方程式を考えない方法 t = 0 1 1 (u[i] - u[i-1])/h 温度u l:熱伝導率 A:断面積 r:密度 c:比熱 h 分割する h h w 1 u[i-1] u[i] u[i+1] (u[i] - u[i-1])/h 温度勾配？ w[i] A B C D E ・・・・・ V 1 時刻 x0 x1 x3 x4 x20 2 3 0.001

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

微分方程式を考えない方法 t = 0 1 1 温度u l:熱伝導率 A:断面積 r:密度 c:比熱 k:時間刻み h h w u[i-1] 分割する h h w 1 u[i-1] u[i] u[i+1] w[i] Q = - (u[i] - u[i-1])/h*lAk Q’ = - (u[i+1] - u[i])/h*lAk w[i] = u[i] +(Q-Q’)/(hArc) = u[i] +(u[i+1]-2u[i]+u[i-1])lk/rch2

エクセルでトライ！ A B C D E ・・・・・ V 1 時刻 x0 x1 x3 x4 x20 2 3 0.001 4 0.002 5 3 0.001 4 0.002 5 0.003 6 時間刻み B2,C2,D3 から導出

エクセルでトライ！ w[i] = u[i] +(u[i+1]-2u[i]+u[i-1])lk/rch2 A B C D E ・・・・・ V 1 時刻 x0 x1 x3 x4 x20 2 3 0.001 4 0.002 5 0.003 6 時間刻み B2,C2,D2 から導出 w[i] = u[i] +(u[i+1]-2u[i]+u[i-1])lk/rch2 lk/rch2=k*(l/rc)/h2=0.001*1/0.05^2=0.4

今日の課題 演習問題６．１ 次の放物型偏微分方程式 を下記の条件で解け． 　　　　　今日の課題 演習問題６．１　次の放物型偏微分方程式 　　　　　　　　　　を下記の条件で解け． ∂u ∂2u 　　　　　= ∂t ∂x2 次回Ｐ７４ 初期条件：t=0, 0<x<1で，u=1 境界条件：t>0, x=0, 1で，u=0 ただし，刻みはh = 1/4, k = 1/32を用いよ．