Хугацааны интеграци /time integration/

Хугацааны интеграци

Тухайн болон бүрэн дифференциал тэгшитгэл хугацаанаас хамаарсан байвал ихэнхдээ хөдөлгөөнтэй холбогдоно. Тухайлбар Нюьтоны 2-р хууль бол уг хөдөлгөөний дифференциал тэгшитгэлийн жишээ юм. Дифференциал тэгшитгэлийг тухайн /partial/ ба бүрэн /ordinary-заримдаа ердийн гэж нэрлэнэ/ гэж ангилахаас гадна шугаман ба шугаман бус, эрэмбээр нь нэг ба хоёр, өндөр эрэмбийн гэх мэт ангилна. Практикт, нэг ба хоёрдугаар эрэмбийн шугаман /linear/ ба шугаман бус /non-linear/ дифференциал тэгшитгэлүүд илүү тохиолдоно. Хугацаанаас хамаарсан дифференциал тэгшитгэлийг шийдэх аналитик аргаас гадна тооцон бодох арга гэж бий. Үүний нэг хялбар жишээ нь төгсгөлөг ялгавараар /Finite difference/ интеграцилах арга юм. Төгсгөлөг ялгаварын ухрах /backward/, давших /forward/, төвийн /central/ гэсэн гурван төрлийн ялгавар байх ба эдгээрээс давших болон төвийн аргыг нь түлхүү ашиглана. Төгсгөлөг ялгаварын аргын тухай, тэдний жишээг урьд нийтлэлүүдэд тусгасан буй.

Товчоор Эйлэрийн, Рунга-Куттагын, Онож заслах буюу Эйлэр-Кошийн гэсэн аргуудыг авч үзье. Жишээнд авсан бодлого:

                                                          

                                                         y'= - 2.3y,   y(0)=1

Эйлэрийн арга

Текстийг дараа оруулъя. Одоохондоо онолыг нь эндээс here уншина уу. Тооцон бодох фортран кодыг сонирхвол:

1:     Program Eiler integration  
2:     integer i  
3:     real*8, dimension (10) :: y1, y2  
4:     real dt  
5:     open(1,file='eiler.dat')  ! output file  
6:     write(*,*) 'please type time step dt'  
7:     read(*,*) dt  
8:     y1(1)=1           ! initial condition  
9:     do i=1,10  
10:      y2(i)=y1(i)-2.3*y1(i)*dt  
11:      y1(i+1)=y2(i)  
12:     write(1,*) y1(i) !,i=1,10)  
13:     end do  
14:     End program  

Бодолтонд хэд хэдэн хугацааны алхамыг 0.1 - ээс эхлэн 0.6 хүртэл авсан болно. Үр дүнгүүдийг аналитик шийдтэй харьцуулж харвал:

Эйлэрийн аргын хувьд хугацааны алхам том байх тусам тохирол байхгүй, тогтворгүй шийдэнд хүргэж байгааг d=0.3-0.6 хоорондох үр дүнгээс харж болно. Харин dt=0.1 байхад бодит шийдтэйгээр дөхөж очиж байгаа боловч яг нарийн тохирохгүй байна. Иймд Эйлэрийн арга нь тохирлын зэрэг /order of accuracy/ муу арга юм. 

Рунга-Куттагийн арга

Энэ аргыг 1900 оны үед Рунга болон М.В.Кутта нар боловсруулсан ба тохирол өндөртэй ил ба далд хосолмол арга юм. Мөн л онолыг одоохондоо эндээс here уншина уу. Тооцон бодох фортран кодыг сонирхвол:

1:     Program Runga Kutta integration  
2:     integer i  
3:     real*8, dimension (10) :: y1, y2   
4:     real dt, k1, k2, k3, k4  
5:     open(1,file='runga.dat')  
6:     write(*,*) 'please type time step dt'  
7:     read(*,*) dt  
8:     y1(1)=1  
9:     do i=1,10  
10:      k1=-2.3*y1(i)  
11:      k2=-2.3*(y1(i)+dt/2*k1)  
12:      k3=-2.3*(y1(i)+dt/2*k2)  
13:      k4=-2.3*(y1(i)+dt*k3)  
14:      y2(i)=y1(i)+(dt/6)*(k1+2*k2+2*k3+k4)  
15:      y1(i+1)=y2(i)  
16:     write(*,*) k2  
17:     write(*,*) y1(i) !,i=1,10)  
18:     write(1,*) y1(i) !,i=1,10)  
19:     end do  
20:     End program  

Энэ кодонд бичигдэж байгаа к - нууд бол Рунга-куттагын аргын өсөлтийн коэффициентууд юм. Хугацааны алхамыг дээрхтэй ижил байдлаар авч бодит шийд болон эйлэрийн хамгийн нарийвчлалтай шийдүүдтэй харьцуулбал:

Зурагт бодит шийд болон РК4 /Рунга-куттагын аргын товчилсон нэр/ -ын хамгийн бага алхамтай үеийн шийдүүд давхцаж харагдаж байгаа боловч тодорхой зөрөө байгаа гэдгийг хэлье. Хэдийн энэ арга өндөр эрэмбийн тохиролтой боловч хугацааны алхамыг буруу сонговол мөн л ойролцоо шийдийг өгч зарим тохиолдолд шийдэл тогтворгүй болох нь бий. 

Эйлэр-кошийн арга

Эйлэр-кошийн арга гэж орос сурах бичигт нэрлэгдэх боловч англи хэл дээрх сурах бичигт бараг л энэ нэрээрээ олдохгүй. Ихэвчлэн Онож заслах /predictor-corrector method/ арга гэдгээр олдоно. Онож заслах аргын тухай монгол ба англи онолыг тус бүр дээр нь дараад орж үзээрэй. Орос номонд дурьдагдах эйлэр-кошийн арга нь жинхэнэ онож заслах аргаас үйлдлийн хувьд 1-ээр дутуу байх боловч агуулга нь маш ойролцоо байдаг. Эхлээд Эйлэр-кошийн гэгдэх онож заслах аргын фортран кодийг авч үзье.

1:     Program Predictor_Corrector integration  
2:     integer i  
3:     real*8, dimension (10) :: y1, y2  
4:     real dt, k1  
5:     open(1,file='pred.dat')  
6:     write(*,*) 'please type time step dt'  
7:     read(*,*) dt  
8:     y1(1)=1  
9:     do i=1,10  
10:      k1=y1(i)+dt*(-2.3*y1(i))  ! onoh alham  
11:      y2(i)=y1(i)+(dt/2)  
12:     1  *(-2.3*y1(i)+(-2.3*k1))  ! zaslah alham  
13:      y1(i+1)=y2(i)  
14:     write(*,*) k1  
15:     write(*,*) y1(i) !,i=1,10)  
16:     write(1,*) y1(i) !,i=1,10)  
17:     end do  
18:     End program  

Энд давхар нарийвчлалтай бодит тоогоор зарласан байгааг анзаарч байгаа байх. Хугацааны алхам мөн дээрхтэй ижил, кодонд дурьдагдсан к нь онох алхам дахь хувьсагчийн утга юм. Ингээд үр дүнг харьцуулж харвал:

Энэ жишээний хувьд РК4 аргатай харьцуулбал илүү тогтвортой боловч тохиролын хувьд бараг л ялгаагүй гэдэг нь ажиглагдсан. Гэхдээ өөр жишээний хувьд авч үзвэл харилцан адилгүй байх нь ойлгомжтой. Ингээд жинхэнэ онож заслах аргын фортран кодыг харвал:

1:     Program Predictor_Corrector integration  
2:     integer i  
3:     real*8, dimension (10) :: y1, y2  
4:     real dt, k1, k2  
5:     open(1,file='pred1.dat')  
6:     write(*,*) 'please type time step dt'  
7:     read(*,*) dt  
8:     y1(1)=1  
9:     do i=1,10  
10:      k1=y1(i)+0.5*dt*(-2.3*y1(i))  !Ehnii onolt  
11:      k2=y1(i)+0.5*dt*(-2.3*k1)    !daraachiin onolt  
12:      y2(i)=2*k2-y1(i)                   !zasah alham  
13:      y1(i+1)=y2(i)  
14:     write(*,*) y1(i) !,i=1,10)  
15:     write(1,*) y1(i) !,i=1,10)  
16:     end do  
17:     End program  

Энэ онож заслах арга нь Эйлэр-кошийн аргыг бодвол хугацааны хагас алхамд хоёр удаа утгыг ондогоороо ялгаатай юм. Эйлэр-кошийн аргын онолыг сүүлд оруулах болно. Ингээд хоёр аргын хувьд ижил хугацааны алхамтай үед хэр зэрэг зөрүүтэй үр дүнг өгдөгийг сонирхоё.

Эхний гэдэг нь Эйлэр-коши, дараа гэдэг нь онож заслах аргын кодоор бодсон үргүүд юм. Хичнээн ойртуулсан ч яг давхацсан мэт график харагдаж байгаа учир зөрүүг тоон байдлаар харъя. 

Ялгаа
dt=0.1/ehnii-daraa/	dt=0.2/ehnii-daraa/
1.49541E-08	-6.64294E-09
5.60705E-08	-1.19385E-08
1.66554E-08	-2.20931E-08
-4.60162E-09	3.98717E-09
1.10899E-08	-7.27205E-09
3.16795E-09	-1.5079E-08
-1.0955E-11	-8.89163E-09
4.59258E-09	-3.39542E-09
3.72246E-09	-1.40696E-09

Маш бага зөрүүтэй байгааг дээрх хүснэгтээс харж болно. Тооны стандарт тэмдэглэгээг уншихад асуудал байвал ийшээ яваарай.

Нэг гадаад гар сонирхсон учир онолгүйгээр шууд кодуудыг орууллаа. Дараа нь нөхөж сайжруулах болно.