Home > Jejak langkah untuk bersekolah lagi, Python, Semester 2, Tutorial > Runge Kutta implisit dengan python (2)

Runge Kutta implisit dengan python (2)

Setelah pada tulisan sebelumnya, telah ada tiga class pendukung yang telah didefinisikan, kita siap untuk mendefinisikan penyelesaian runge kutta implisit. Skenarionya, digunakan dua metode iterasi, masing-masing Newton dan Picard. Metode Newton akan dibahas lebih dulu.
Berdasarkan dua persamaan yang dijelaskan pada tulisan pertama, serta metode penyelesaian yang dijelaskan dalam [1], persamaan nilai antara diolah sedemikian rupa sehingga terbentuk persamaan berikut. J adalah fungsi Jacobi untuk setiap persamaan yang terlibat.

J \delta = -F(y)      (1)

Selanjutnya, \delta diselesaikan sehingga nilainya lebih kecil atau sama dengan nilai toleransi. Saat kondisi \delta tercapai, nilai \delta digunakan untuk memperbaiki nilai antara menggunakan persamaan (2). Yang perlu diperhatikan adalah elemen \delta harus sebanyak jumlah persamaan dikalikan orde penyelesaian.

y_{j+1} = y_{j} + \delta      (2)

Nilai antara yang diperoleh dari persamaan di atas, digunakan untuk memperperbaiki nilai fungsi di titik step size dikalikan k. k adalah nilai step ke sekian dalam iterasi. Nilai fungsi diperbaiki dengan persamaan (3).

x_{k+1}=x_{k} + h \Sigma_{j=1}^{s} b_{i}f(t_{k}+c_{j}h,y_{ik}      (3)

Nilai fungsi yang baru digunakan kembali untuk memprediksi nilai antara, melalui persamaan (4).

y_{ik}=x_{k}+h \Sigma_{j=1}^{s} a_{ij}f(t_{k}+c_{j}h,y_{jk})      (4)

Jika dikaitkan dengan persamaan (1), nilai F(y) diperoleh dari persamaan (4), tepatnya seperti persamaan (5).

F(y) = y_{ik} - x_{k} - h \Sigma_{j=1}^{s} a_{ij}f(t_{k}+c_{j}h,y_{jk})      (5)

Sedangkan nilai J dari persamaan (1), merupakan matriks jacobian dari persamaan yang terlibat dalam sistem. Sebagai contoh, sistem dengan dua persamaan dan akan diselesaikan dengan RK orde dua, maka akan terdapat matriks J berukuran 4 x 4. Kenapa bisa begitu?

Dari faktor orde penyelesaian, akan diperoleh matriks berdimensi orde x orde. Nah, setiap elemen dalam matriks berdimensi orde x orde tersebut, sejatinya merupakan matriks juga yang berdimensi n x n, dengan n adalah banyaknya persamaan (dan sudah tentu banyaknya variabel) yang terlibat dalam sistem. Karena orde penyelesaian yang digunakan adalah dua, dan jumlah persamaannya juga dua, maka dimensi matriks keseluruhan adalah 4 x 4. Kesimpulannya, sistem dengan m persamaan dan diselesaikan dengan orde n, akan dihasilkan matriks J dengan dimensi mn x mn. Lalu bagaiman mengisi setiap elemen matriks tersebut?

Matriks pertama berdimensi orde x orde, telah disebutkan berisi matriks berdimensi n x n, dengan n adalah jumlah persamaan yang terlibat. Untuk matriks yang kedua ini, isinya adalah sebagai berikut. Jika terdapat sistem dengan dua persamaan, masing-masing y_{1} dan y_{2}, serta variabel masing-masing var_{1} dan var_{2}, maka isi matriks kedua adalah seperti persamaan (6)

\begin{bmatrix}  \frac{\partial{y_{1}}}{\partial{var_{1}}} & \frac{\partial{y_{1}}}{\partial{var_{2}}}\\  \frac{\partial{y_{2}}}{\partial{var_{1}}} &  \frac{\partial{y_{2}}}{\partial{var_{2}}}  \end{bmatrix}       (6)

Jika orde penyelesaian sama dengan dua, maka persamaan (6) adalah bagian dari elemen ke i,j dari matriks berdimensi orde x orde. Matriks berdimensi 2 x 2 tersebut adalah seperti persamaan (7), dengan I adalah matriks identitas, sedangkan j_{i,k} adalah matriks dari persamaan (6)

\begin{bmatrix}  I-ha_{1,1}j_{1,1} & -ha_{1,2}j_{1,2} \\  ha_{2,1}j_{2,1} & I-ha_{2,2}j_{2,2}  \end{bmatrix}       (7)

Terkahir, program python untuk menyelesaikan sistem persamaan non-linier dengan runge kutta implisit adalah seperti program berikut ini. Program tersebut mengunakan asumsi bahwa ketika dalam suatu step, iterasi untuk mencapai konvergensi (nilai \delta lebih kecil daripada nilai toleransi) melebihi 5000, maka iterasi dihentikan. Perhitungan runge kutta dengan iterasi Newton selanjutnya dihentikan.


from numpy import *
from numpy import linalg as la
import os.path
from time import *
from Equations import *
from Butcher import *
from Symbols import *
class NewtonRKSolver():
	def __init__(self,Eq,Symbols,Butcher,Initial,Stop,StepSize,TimeTarget):
		"""lenEq: banyaknya persamaan"""
		self.lenEq=len(Eq)
		
		"""persamaan yang akan diselesaikan"""
		self.equations=Eq
		
		"""simbol/variabel yang digunakan"""
		self.symbols=Symbols
		
		"""orde: orde penyelesaian"""
		self.orde=len(Butcher[1])
		
		"""Jv=-F: sistem persamaan linier yang akan diselesaikan"""
		self.F=arange(1)
		self.v=arange(1)
		self.J=arange(1)
		
		"""x: nilai target yang akan dihitung pada setiap step size"""
		self.x1=array(Initial)
		
		self.errorStatus=0
		
		"""
		self.stop: stopping criteria
		self.h=step size
		"""
		self.stop=Stop
		self.h=StepSize
		self.t=TimeTarget
		self.selisih=self.t/self.h-int(self.t/self.h)

		"""y: nilai antara dalam bentuk sub array untuk menghitung derivatif"""
		"""y1: nilai antara, yang tebakan awalnya disamakan dengan nilai awal"""
		tempY=[]
		for i in range(self.orde):
			tempY.append(Initial)
		self.y=array(tempY)	
		self.y1=array(tempY)
		self.y1=reshape(self.y1,(self.orde*self.lenEq,1))
		l=0

		"""delta: perbaikan nilai antara"""
		self.delta=self.y1
		
		"""a,b,c: komponen array Butcher"""
		self.a=Butcher[0]
		self.b=Butcher[1]
		self.c=Butcher[2]
		
		"""membentuk vektor F"""
		self.__setF()


		"""membentuk matriks J"""
		self.__setJ()
		
		self.elapsed=0
		self.start=time()
		if self.selisih==float(0):
			i=0
			while(i<int(self.t/self.h) and self.errorStatus==0):
				self.__iterate()
				self.__setY()
				self.__setF()
				self.__setJ()
				i+=1
			self.elapsed=time()-self.start	

		else:
			i=0
			while(i<int(self.t/self.h) and self.errorStatus==0):
				self.__iterate()
				self.__setY()
				self.__setF()
				self.__setJ()
				i+=1
			self.h=self.selisih
			self.__iterate()
			i=int(self.t/self.h)+1
			self.elapsed=time()-self.start

		if self.errorStatus==1:
			print 'Quit unseccessful in ' + str(self.elapsed) + 'seconds ... '
			print 'Tidak konvergen di tahap-' + str(i)
		
	def __getNumDerivative(self,eq,number,val):
		a=self.equations[eq]
		b=self.equations[eq]
		for i in range(self.lenEq):
			if i==number:
				a=a.subs(self.symbols[i],(val[i]+0.001))
				b=b.subs(self.symbols[i],val[i])
			else:
				a=a.subs(self.symbols[i],val[i])
				b=b.subs(self.symbols[i],val[i])
		return (a-b)/0.001

	def __setF(self):
		p=[]
		l=0
		for i in range(self.orde):
			for j in range(self.lenEq):
				q=0
				for k in range(self.orde):
					q=q+self.a[i][k]*self.y[k][j]
				p.append(self.x1[j]-self.y1[l]+q*self.h)
				l+=1
		self.F=array(p)
		
		self.F=reshape(self.F,(self.orde*self.lenEq,1))
		
	def __setJ(self):
		tempJ=[]	
		for i in range(self.orde):
			for l in range(self.orde):
				tempJ1=[]
				for j in range(self.lenEq):
					for k in range(self.lenEq):
						tempJ1.append(self.__getNumDerivative(j,k,self.y[l]))
				tempJ1=array(tempJ1)
				tempJ1=reshape(tempJ1,(self.lenEq,self.lenEq))
				tempJ.append(tempJ1)
		tempJ=array(tempJ)
		p=[]
		k=0
		for i in range(self.orde):
			for j in range(self.orde):
				if i==j:
					q=identity(self.lenEq)-self.h*self.a[i][j]*tempJ[k]
					p.append(q)
				else:
					q=-self.h*self.a[i][j]*tempJ[k]
					p.append(q)
				k+=1
		self.J=array(p)
		p=[]
		for i in range(self.orde):
			for j in range(self.lenEq):
				for k in range(self.orde*i,self.orde*(i+1)):
					for l in range(self.lenEq):
						p.append(self.J[k][j][l])
		self.J=array(p)
		x=self.orde*self.lenEq
		self.J=reshape(self.J,(x,x))
	
	def __setX(self):
		self.y1=reshape(self.y1,(self.orde,self.lenEq))
		for i in range(self.lenEq):
			tempY=0
			for j in range(self.orde):
				tempY=tempY+(self.b[j]*self.y1[j][i])
			self.x1[i]=self.x1[i]+(self.h*tempY)
		self.y1=reshape(self.y1,(self.orde*self.lenEq,1))
		
	def __setY(self):
		l=0
		for i in range(self.orde):
			for j in range(self.lenEq):
				self.y[i][j]=self.y1[l]
				l+=1
				
	def __iterate(self):
		l=0
		self.delta=linalg.solve(self.J,self.F)
		self.y1=self.delta+self.y1
		l+=1
		while(self.stop<linalg.norm(self.delta) and l=5000:
			self.errorStatus=1
			
	def getF(self):
		return self.F
		
	def getJ(self):
		return self.J

Referensi
[1]. G. Hall and J.M. Watt. Modern Numerical Methods for Ordinary Differential Equations. Oxford University Press, 1976.

About these ads

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

%d bloggers like this: