Studi Aplikasi Metode Operator Splitting untuk Menyelesaikan Anisotropic Problem
Is Bunyamin Suryo a, *, Maureen Clerc b
a Departemen Teknik Mesin, Fakultas Teknologi Industri, ITS
ITS Kampus Keputih Sukolilo Surabaya 60111
b Athena Project Team - INRIA INRIA Sophia Antipolis, Prancis *E-mail: isbunyamin@gmail.com

## Abstract
Electroencephalography  (EEG)  is  a  non-invasive  technique  for  studying  electrical  activity  in  the electrical  brain.  The  electrical  activity  is  a  complex  process  of  electrical  propigation  due  to  the complex structure of the brain. This complex brain structure causes parts of the brain to have different conductivity, both in size and orientation, which is called anisotropic conductivity. Using Maxwell's Equation  approach,  electrical  activity  in  the  brain has been  intensively  studied.  For  simplicity,  the quasistatic Maxwell Equation (Poisson's Equation) is used to model electrical activity in the brain. In this study, an initial study was conducted using a new method, Operator Splitting Method (OSM), to solve  anisotropic  (2D)  2-Dimensional  Poisson  Equation.  A  freeware  of  the  finite  element  method (FEM) is used to construct the matrix used in the OSM algorithm. The OSM algorithm is then tested to solve anisotropic Laplace 2D Equation and anisotropic Poisson Equation with a dipolar source. After that, the OSM solution is validated using theoretical solutions and direct numerical solutions. Using L2-Error Norm, the convergence rate of the OSM algorithm is then analyzed. Several numerical experiments have been carried out to test the OSM algorithm. The OSM solution of the anisotropic Laplace 2D Equation corresponds to a theoretical solution and a direct numerical solution.  For 2D anisotropic Poisson Equations with dipolar sources, several similar results have been obtained. The OSM  solution  pattern  is  very  close  to  the  numerical  solution  pattern  directly.  The  results  of  this numerical experiment raise hopes to try to apply the OSM algorithm to more complex problems such as real modeling on the human head.
Keywords: Maxwell's Equation, FEM, Laplace's Equation, Poisson's Equation

## Abstrak
Electroencephalography (EEG) adalah  teknik  non-invasif  untuk  mempelajari  aktivitas  listrik  pada otak listrik. Aktivitas listrik tersebut adalah proses propigasi listrik yang komplek karena struktur otak yang rumit. Struktur otak yang komplek ini menyebabkan bagian-bagian otak memiliki konduktivitas yang  berbeda,  baik  besarnya  maupun  orientasinya,  yang  disebut  konduktivitas  anisotropik. Menggunakan  pendekatan  Persamaan  Maxwell,  aktivitas  listrik  pada  otak  telah  dipelajari  secara intensif.  Untuk  penyederhanaan,  Persamaan  Maxwell  quasistatik  (Persamaan  Poisson)  digunakan untuk memodelkan aktivitas listrik pada otak. Dalam penelitian ini, dilakukan studi awal penggunaan metode baru, Operator Splitting Method (OSM), untuk menyelesaikan Persamaan Poisson 2-Dimensi anisotropik (2D). Sebuah freeware dari metode elemen hingga (FEM) digunakan untuk membangun matriks yang digunakan dalam algoritma OSM. Algoritma OSM kemudian diuji untuk menyelesaikan Persamaan 2D Laplace anisotropik dan Persamaan Poisson anisotropik dengan sumber dipolar. Setelah itu, solusi OSM divalidasi dengan menggunakan solusi teoritis dan solusi numerik langsung. Dengan menggunakan L2-Error Norm , laju konvergensi dari algoritma OSM kemudian dianalisis. Beberapa percobaan numerik telah dilakukan untuk menguji algoritma OSM. Solusi OSM dari Persamaan 2D Laplace anisotropik bersesuaian dengan solusi teoritis dan solusi numerik langsung. Untuk Persamaan 2D Poisson anisotropik dengan sumber dipolar, beberapa hasil serupa juga telah diperoleh. Pola solusi OSM  sangat  mendekati  pola  solusi  numerik  langsung.  Hasil  percobaan  numerik  ini  memunculkan harapan  untuk  mencoba  menerapkan  algoritma  OSM  untuk  masalah  yang  lebih  kompleks  seperti pemodelan riil pada kepala manusia.
Kata kunci: Persamaan Maxwell, FEM, Persamaan Laplace, Persamaan Poisson

## 1. Pendahuluan
Electroencephalography (EEG)  adalah  teknik  non-invasif  untuk  mempelajari  aktivitas  listrik  pada  otak  listrik. Aktivitas listrik tersebut adalah proses propigasi listrik yang komplek karena struktur otak yang rumit. Struktur otak yang komplek ini menyebabkan bagian-bagian otak memiliki konduktivitas yang berbeda, baik besarnya maupun orientasinya, yang disebut konduktivitas anisotropik. Menggunakan pendekatan Persamaan Maxwell, aktivitas listrik pada otak telah dipelajari secara intensif. Untuk penyederhanaan, Persamaan Maxwell quasistatik (Persamaan Poisson) digunakan untuk memodelkan  aktivitas  listrik  pada  otak,  seperti  ditunjukkan  oleh  Persamaan  1.  Perdekatan  ini  sangat  valid  untuk pemodelan EEG karena frekuensi EEG pada kepala manusia sangat lemah, lebih rendah dari 10 2 Hz [1].
  p J V        (1)
dimana V adalah potensial listrik, J p adalah arus utama, dan σ adalah sifat konduktivitas. Dalam pemodelan EEG pada kepala manusia, sifat konduktivitas Persamaan 1 harus ditentukan. Properti ini terkait dengan struktur kepala manusia. Struktur  kepala  manusia  adalah  struktur  yang  sangat  kompleks  karena  merupakan  struktur  heterogen  dengan  lapisan jaringan yang berbeda dan setiap lapisan jaringan memiliki sifat khusus seperti yang ditunjukkan oleh Gambar 1 [2-3]. Untuk  pemodelan  kepala  manusia  yang  realistis,  sangat  sulit  untuk  menyelesaikan  model  yang  mempertimbangkan kompleksitas struktur kepala manusia [4-6]. Secara umum, model ini disebut model kepala manusia anisotropik. Ide dari penelitian  ini  adalah  untuk  mengusulkan  metode  baru,  yang  disebut  Operator  Splitting  Method  (OSM),  untuk memecahkan masalah anisotropik. Metode yang diusulkan diharapkan dapat menyelesaikan masalah anisotropik dengan mudah.  Dalam  penelitian ini,  studi  kelayakan menggunakan  OSM  untuk memecahkan masalah  anisotropik  sederhana akan disajikan. Harapannya adalah studi kelayakan ini dapat memotivasi para peneliti untuk menerapkan OSM dalam model kepala manusia yang realistis.
Banyak peneliti telah mencoba mengembangkan model kepala manusia untuk menyelesaikan forward problem of EEG source localization . Forward problem of EEG source localization bertujuan untuk menghitung medan listrik yang dihasilkan  oleh  arus  primer  pada  dalam  geometri/domain  tertentu.  Metode  elemen  batas  (BEM)  dan  metode  elemen hingga  (FEM)  dapat  diimplementasikan  untuk  menyelesaikan forward  problem  of  EEG  source  localization .  Kedua metode  pendekatan  numerik  ini  memiliki  kelebihan  dan  kekurangan.  Simulasi  dengan  pendekatan  BEM  memerlukan biaya komputasi yang relative rendah karena dekomposisi domain atau meshing hanya diterapkan pada batas domain, seperti yang ditunjukkan pada Gambar 2a.
(a)                                                   (b) Gambar 1. (a) The MRI of the human head (b) The human head tissues [7] .
Kelemahan dari pendekatan BEM adalah pada asumsi yang digunakan sehingga menjadi tidak realistis, yaitu semua lapisan jaringan otak dianggap homogen dan isotropik [8]. Sedangkan pendekatan dengan simulasi FEM membutuhkan biaya  komputasi  yang  tinggi  karena  seluruh  domain  komputasi  harus  didekomposisi,  seperti  yang  ditunjukkan  oleh Gambar  2b.  Ini  berarti  bahwa  jumlah  node  atau  titik  yang  tidak  diketahui  pada  domain  jauh  lebih  banyak  daripada menggunakan BEM. Namun, pendekatan FEM jauh lebih realistis daripada pendekatan BEM karena kompleksitas lapisan jaringan  otak  yang heterogen dan anisotropik  dapat  diakomodasi.  Dalam  penelitian  ini,  menggunakan  OSM,  masalah anisotropik  sederhana  dicoba  untuk  dipecahkan  dengan mempertimbangkan model  isotropiknya  dan memperkenalkan faktor koreksi, yang disebut lineic.
(a)                                                                          (b) Gambar 2. (a) Human head model through triangle BE (b) Human head model through tetrahedral FE .

## 2. Metode Penelitian
Pertama,  OSM  akan digunakan  untuk menyelesaikan  Persamaan  2D  Laplace anisotropik.  Solusi  OSM  kemudian dibandingkan dengan solusi numerik langsung dan solusi analitis untuk memvalidasi solusi OSM. Setelah itu, OSM akan diuji untuk mendapatkan solusi Persamaan 2D Poisson anisotropik. Berdasarkan Persamaan Poisson (Persamaan 1), bila sisi kanan Persamaan 1 sama dengan nol, maka Persamaan 1 merupakan Persamaan Laplace, seperti yang dinyatakan oleh  Persamaan  2.  Pada  domain  2D,  operator  konduktifitas ( delta )  pada  Persamaan  1  dan  2 dapat  dinyatakan  sebagai Persamaan 3.
  0     V  (2)
       y x    0 0 (3)
dimana x  adalah konduktifitas arah-x, y  adalah konduktifitas arah-y. Bila y x    , maka  adalah konduktifitas isotropik, iso  . Bila 0   y x   atau y x     0 , maka  adalah konduktifitas anisotropik, anis  . Sehingga, 2D Poisson’s and Laplace’s equation dapat dinyatakan secara berturut-turut seperti Persamaan 4 dan 5.
  p y x J y V x V V             2 2 2 2    (4)   0 2 2 2 2           y V x V V y x    (5)
Untuk mengimplementasikan OSM, operator konduktifitas pada anisotropic 2D Poisson’s - Laplace’s equation dianggap sebagai penjumlahan dari isotropic conductivity , iso  , dan lineic/correction conductivity , lin  , seperti berikut:
lin iso aniso      (6)
Bila diasumsikan y x    , Persamaan 6 dapat ditulis sebagai berikut:
         
              0 0 0 0 0 0 0 lin x y y lin iso y x anis         (7)
Bila Persamaan 6 disubtitusikan ke Persamaan 1 dan 2, dengan asumsi 1  y  , Persamaan 1 dan 2 dapat disusun ulang secara berturut-turut sebagai Persamaan 8 dan 9. Nilai dari lin x   tergantung pada seberapa kuat x  terhadap y  (rasio antara konduktifitas ke arah-x dan arah-y).
p lin x y x J x V y V x V y V x V                         2 2 2 2 2 2 2 2 2 2    (8) 0 2 2 2 2 2 2 2 2 2 2                       x V y V x V y V x V lin x y x   
(9)
Menyelesaikan  Persamaan  Poisson dengan  pendekatan  FEM  akan menghasilkan  sebuah  sistem linier B AV  . Pada dasarnya, sistem linier B AV  dibangun berdasarkan proses sebagai berikut: matrik A dibangun dari left hand- side of the equation , sedangkan matrik B dibangun dari right hand-side of the equation . Penggunaan OSM, ada 2 term di left hand-side of the equation . Matrik pertama A dibangun dari            2 2 2 2 y V x V , disebut sebagai iso A . Matrik kedua
A dibangun dari 2 2 x V lin x     , disebut sebagai lin A . Dengan freeware FreeFEM++ [9], kedua term tersebut didiskritisasi sehingga  dihasilkan  elemen  matrik  A.  Setelah  itu,  algoritma  OSM  akan  dikembangkan  berdasarkan  matrik  A  yang dihasilkan sebelumnya. Pada penelitian ini, algoritma OSM yang ditawarkan adalah sebagai berikut:
A iso .V n = B - A lin .V (n-1) (10)
V (n+1) =  .V n + (1-  ).V (n-1) (11)
dimana  adalah koefisien relaksasi dan n adalah jumlah iterasi. Koefisien relaksasi ini digunakan untuk mempercepat konvergensi  saat  proses  iterasi.  Algoritma  OSM  yang  ditulis  dengan  Matlab  kemudian  diuji  untuk  menyelesaikan anisotropic Laplace’s - Poisson’s equation with dipolar source . Dengan menggunakan L2-Error Norm , laju konvergensi algoritma OSM akan dianalisa.

## 3. Hasil dan Pembahasan

## 3.1 Solusi Numerik dari Persamaan 2D Laplace Anisotropik
Sebelum  menyelesaikan  Persamaan  2D  Poisson  anisotropik,  OSM  akan  diuji terlebih  dulu  untuk  menyelesaikan Persamaan 2D Laplace anisotropik (Persamaan 12) dimana
y x   3  .   0 2 2 2 2           y V x V V y x    (12)
Domain  dari  simulasi  ini  adalah  sebuah  persegi  panjang  pada  koordinat-xy  seperti  ditunjukan  oleh  Gambar  3.  Batas kondisi untuk menyelesaikan Persamaan 12 adalah sebagai berikut:
  0 ,  y a V (13)   0 0 ,    x n V ,   0 ,    b x n V ,   1 , 0     y n V
(14)
Gambar 4. Domain decomposition/meshing .
Pada  batas  domain,  ada  dua  jenis  kondisi  batas  yang  digunakan,  yaitu  batas  kondisi  Dirichlet  dan  batas  kondisi Neumann. Sebuah freeware , FreeFEM ++, digunakan untuk mendiskritisasi domain simulasi dan membangun matrik dari sistem linier B AV  .  Kemudian algoritma  OSM  yang  ditulis  dengan  Matlab akan  digunakan  untuk menyelesaikan sistem linear ini dan menampilkan solusi numerik. Dimensi domain persegi panjang 2D adalah (0,1x0,7), didekomposisi menggunakan meshing segitiga menjadi 1852 node atau simpul dan 3462 segitiga seperti yang ditunjukkan oleh Gambar 4. Solusi analitis atau teoritis dan solusi dari OSM dari Persamaan 12 ditunjukkan oleh Gambar 5 dan 6 secara berturut- turut. Hasil analitis diperoleh dengan menggunakan metode Pemisahan Variabel, sedangkan solusi dari OSM diperoleh dengan menggunakan koefisien relaksasi (Persamaan 11) 0,1 dengan 30 iterasi. Gambar 7 menunjukkan perbandingan antara solusi analitis dengan OSM. Dari Gambar 5, 6, dan 7 dapat dilihat bahwa solusi yang dihasilkan OSM memiliki kesamaan ( good agreement ) dengan hasil analitis. Laju konvergensi algoritma OSM dapat dilihat pada Gambar 8. Secara umum,  algoritma  yang  diusulkan  bekerja  dengan  baik  karena  tren  konvergensi  menurun  sepanjang  iterasi.  Dalam penelitiannya, Suryo [10] menjelaskan pengaruh dari nilai koefisien relaksasi pada kecepatan konvergensi dari simulasi OSM.  Suryo  [10]  menvariasikan  nilai  koefisien  relaksasi  dari  0,03  hingga  0,6.  Dengan  nilai  koefisien  relaksasi  yang semakin besar maka konvergensi semakin cepat dicapai. Dari hasil ini, OSM bisa diuji untuk menyelesaikan persamaan yang lebih komplek, yaitu Persamaan 2D Poisson anisotropik dengan sumber dipolar.
Gambar 5. Analytical solution of anisotropic 2D Laplace’s equation for
y x   3  .
Gambar 6. OSM solution of anisotropic 2D Laplace’s equation for
y x   3  .
Gambar 7. Comparison between exact and OSM solution in solving anisotropic 2D Laplace’s equation.

## 3.2 Solusi Numerik dari Persamaan 2D Poisson Anisotropik
Pada bagian ini, Persamaan 2D Poisson anisotropik dengan sumber dipolar akan diselesaikan dengan menggunakan OSM.  Untuk  percobaan,  Persamaan  2D  Poisson  anisotropik  dengan  sumber  dipolar  ditentukan  di  mana y x   8  . Solusi ini diperoleh dengan koefisien relaksasi = 0,1 dengan 30 iterasi. Gambar 9 menunjukan grafis hasil simulasi dengan OSM, sedangkan Gambar10 menunjukkan perbandingan antara solusi dengan OSM dan solusi numerik langsung. Dari Gambar 10 dapat disimpulkan bahwa perbandingan antara solusi OSM dengan solusi numerik langsung kesamaan yang cukup baik ( good agreement ).
Gambar 8. Convergence rate of OSM in solving anisotropic 2D Laplace’s equation for
1 . 0   , n = 30.
Gambar 9. Operator Splitting solution of anisotropic 2D Poisson’s equation.
Gambar 10. Comparison between direct numerical and OSM solution in solving anisotropic 2D Poisson’s equation.

## 4. Kesimpulan
Eksperimen numerik  yang  telah  dilakukan  pada penelitian ini  menunjukkan  bahwa  OSM  yang  diajukan mampu menyelesaikan dengan baik Persamaan 2D Laplace anisotropik maupun Poisson anisotropik dengan sumber dipolar. Hasil simulasi OSM memiliki kesamaan yang cukup baik (good agreement) dengan solusi analitis/teoritis maupun dengan hasil solusi numerik langsung. Kinerja algoritma OSM yang diusulkan bekerja dengan baik dalam menyelesaikan Persamaan 2D Laplace anisotropik atau Persamaan 2D Poisson anisotropik dengan sumber dipolar. Metode ini masih perlu banyak pengujian, misalnya dengan menvariasikan besarnya koefisien relaksasi, rasio konduktifitas ke arah-x dank e arah-y. Dari penelitian atau studi awal ini, ada harapan atau motivasi untuk menerapkan algoritma OSM untuk masalah yang lebih kompleks  seperti  model  kepala  manusia  yang  realistis.  Menerapkan  algoritma  OSM  dalam  domain  yang  sedemikian kompleks,  tugas  yang  paling  sulit  adalah  membangun  matriks  yang  akan  diterapkan  pada  algoritma  yang  diusulkan. Dengan demikian di masa yang akan dating, hal ini akan menjadi topik penelitian yang sangat menarik bagi para peneliti. Semoga, metode baru ini (OSM) akan mampu membantu dalam memecahkan masalah pemodelan EEG yang realistis sesuai dengan struktur anatomi otak manusia.

## Daftar Pustaka
[1]  Clerc, M and Papadopoulo, T., 2010, “ Inverse problems in functional brain imaging ”, Reference material for the MVA Master 2 class”.
[2]  Geddes,  L.,  Baker,  L.,  1967,  “The  specific  resistance  of  biological  material.  A  compendium  of  data  for  the biomedical engineer and physiologist”, Medical and biological engineering 5, 271– 293.
[3]  Haueisen,  J., 1996,  “ Methods  of  numerical  field  calculation  for  neuromagnetic  source localization ”, PhD thesis, Shaker-Verlag Aachen, ISBN 3-8265-1691-5.
[4]  Rush, S., Driscoll, D., 1968, “Current distribution in the brain from surface electrodes”, Anesthesia & Analgesia 47 (6), 717– 723.
[5]  Akhtari, M., Bryant, H., Marmelak, A., Flynn, E., Heller, L., Shih, J., Mandelkern, M., Matlachov, A., Ranken, D., Best,  E.,  DiMauro,  M.,  Lee,  R.,  Sutherling,  W.,  2002,  “Conductivities  of  three-layer  live  human  skull”, Brain Topography 14 (3), 151– 167.
[6]  Nicholson, P., 1965, “Specific impedance of cerebral white matter”, Experimental Neurology 13, 386-401.
[7]  Haueisen,  J.,  Tuch  D.S.,  Ramon,  C.,  Schimpf,  P.,  Wedeen,  V.J.,  George,  J.S.  and  Belliveau,  J.W.,  2002,  “The influence of brain anisotropy on human EEG and MEG”, NeuroImage , 15: 159-166.
[8]  Awada, Kassem A., et.al., 1997, “Computational Aspects of Finite Element Modeling in EEG Source Localization”, IEEE Transaction on Biomedical Engineering , Vol. 44, No. 8.
[9]  http://www3.freefem.org/
[10]  Suryo, Is Bunyamin, 2010, “ Operator Splitting for Solving Anisotropic Problem: a Feasibility Study ”, Master thesis, EPU Nice-Sophia Antipolis, Prancis.