Transcript 数値計算法I - 九州工業大学

数値計算法I
第10回
連立微分方程式
2006年度
九州工業大学工学部電気電子工学コース用講義資料
講師：趙孟佑
[email protected]
1
高階微分方程式
• m階の常微分方程式
m
2
m1
d y
dy d y
d y
 f (t, y, , 2 ,L
)
m
m 1
dt
dt dt
dt
これを連立微分方程式に直して、ルンゲ・クッタ法
で解く。
特に2階微分方程式は多くの物理現象で使う。
例:LCR回路
2
高階微分方程式
• m階の常微分方程式
m
2
m1
d y
dy d y
d y
 f (t, y, , 2 ,L
)
m
m 1
dt
dt dt
dt
z1 (t)  y(t), z2 (t)  y(t),L , zm (t)  y(m1) (t) とおくと、
dz1
dz2
dzm1
dzm
(m1)
 y   z2 ,
 y  z3 ,L ,
y
 zm ,
 y(m)
dt
dt
dt
dt
である。
3
連立微分方程式
dmy
dy d 2 y
d m1 y
 f (t, y, , 2 ,L
)
m
m 1
dt
dt dt
dt
を解くことは、下の1階連立微分方程式を解くに等しい
 dz1
 dt  z2

 dz2  z3
 dt

M
 dz
 m1  zm
 dt
 dzm
(m )

y
 f (t, z1 , z2 , z3 ,L , zm )
 dt
4
連立微分方程式
 dy1
 dt  f1 (t, y1 , y2 ,L , ym )

 dy2  f2 (t, y1 , y2 ,L , ym )
 dt

M
 dy
 m1  fm1 (t, y1 , y2 ,L , ym )
 dt
 dym
 dt  fm (t, y1 , y2 ,L , ym )
t=t0での初期条件は
y1 (t0 )  b1, y2 (t0 )  b2 ,L , ym1 (t0 )  bm1, ym (t0 )  bm
5
連立微分方程式の例
交流電界中に置かれた電子の速度：
dv
e

E0 sin  t
dt
me
電子の座標x(t)は
dx
v
dt
よって、t=t0での初期状態をx(t0)=x0,v(t0)=v0として、
 dx
 dt  v
 dv
e
 
E0 sin  t
me
 dt
を解いてx(t),v(t)を求める。
6
2階連立微分方程式のルンゲ・クッタ法による解法
 dy
 dt  f (t, y, z)

 dz  g(t, y, z)
 dt
t  ti で y  yi , z  zi とする。
t=ti+1=ti+hでの値、yi+1,zi+1は次式で与えられる。
h

y

y

k1  2k2  2k3  k4 

i
 i1
6

 zi1  zi  h l1  2l2  2l3  l4 

6
7
2階連立微分方程式のルンゲ・クッタ法による解法
 k1  f (ti , yi , zi )
l  g(t , y , z )
i
i i
1
h
h
h

k

f
(t

,
y

k
,
z

l1 )
i
i
1 i
 2
2
2
2

l2  g(ti  h , yi  h k1 , zi  h l1 )

2
2
2

 k3  f (ti  h , yi  h k2 , zi  h l2 )
2
2
2


h
h
h
l3  g(ti  , yi  k2 , zi  l2 )
2
2
2

 k4  f (ti  h, yi  hk3 , zi  hl3 )
l  g(t  h, y  hk , z  hl )
i
i
3 i
3
4
8
2階常微分方程式の例
2
d x
dx
 3  2x  0
2
dt
dt
初期条件はt=0で、x=0,dx/dt=-1
1階の連立微分方程式に直す。
 dx
 dt  y

 dy  3y  2x
 dt
初期条件はt=0で、x=0,y=-1
刻み幅h=0.2でt=0から1.0まで計算
x(t)  e2t  et
9
2階常微分方程式の例
2
d x
dx
 3  2x  0
2
dt
dt
解析解は x(t)  e2t  et
t
0.00000000
0.20000000
0.40000000
0.60000000
0.80000000
1.00000000
xi(t)
0.00000000
-0.14833333
-0.22088811
-0.24751482
-0.24734177
-0.23246922
解析解
誤差
0.00000000
0.00000000
-0.14841071 0.00007737
-0.22099108
0.00010297
-0.24761742 0.00010260
-0.24743245 0.00009068
-0.23254416 0.00007493
10
2階常微分方程式の例
implicit real*8 (a-h,o-z)
parameter(nmax=100000)
dimension t(nmax),x(nmax),y(nmax)
real*8 k1,k2,k3,k4
real*8 l1,l2,l3,l4
c
c
h=0.2d0
n=1.0d0/h
t(1)=0
x(1)=0
y(1)=-1
do i=1,n
t1=t(i)
x1=x(i)
y1=y(i)
k1=f(t1,x1,y1)
l1=g(t1,x1,y1)
c
t2=t(i)+h/2
x2=x(i)+h/2*k1
y2=y(i)+h/2*l1
k2=f(t2,x2,y2)
l2=g(t2,x2,y2)
c
t3=t(i)+h/2
x3=x(i)+h/2*k2
y3=y(i)+h/2*l2
k3=f(t3,x3,y3)
l3=g(t3,x3,y3)
t4=t(i)+h
x4=x(i)+h*k3
y4=y(i)+h*l3
k4=f(t4,x4,y4)
l4=g(t4,x4,y4)
function f(t,x,y)
implicit real*8 (a-h,o-z)
f=y
return
end
c
function g(t,x,y)
implicit real*8 (a-h,o-z)
g=-2*x-3*y
return
end
c
x(i+1)=x(i)+(k1+2*k2+2*k3+k4)/6*h
y(i+1)=y(i)+(l1+2*l2+2*l3+l4)/6*h
t(i+1)=t(i)+h
end do
do i=1,n+1
sol=dexp(-2*t(i))-dexp(-t(i))
write(6,100)t(i),x(i),sol,dabs(x(i)-sol)
end do
100 format(4(1x,f15.8))
stop
end
c
11