Week-05 (Finite Difference Methods)

Kembali ke Persamaan Diferensial Numerik

yang akan dibahas - Metode Beda Hingga untuk Masalah Linear

Metode Beda Hingga untuk Masalah Nonlinear

Metode ini digunakan untuk mengaproksimasi masalah linear dalam bentuk:

\[\begin{gathered} y^{\prime \prime}=p(x) y^{\prime}+q(x) y+r(x), \quad a \leq x \leq b \\ y(a)=\alpha, y(b)=\beta \end{gathered}\] \[\begin{gathered} w_{0}=\alpha, \quad w_{n+1}=\beta \\ -\left(1+\frac{h}{2} p\left(x_{i}\right)\right) w_{i-1}+\left(2+h^{2} q\left(x_{i}\right)\right) w_{i}-\left(1-\frac{h}{2} p\left(x_{i}\right)\right) w_{i+1}=-h^{2} r\left(x_{i}\right) \end{gathered}\]

Bentuk tersebut dapat dibuat sebagai suatu SPL:

\[ A \mathbf{w}=\mathbf{b} \]

\(\mathbf{w}=\left[\begin{array}{c}w_{1} \\ w_{2} \\ \vdots \\ w_{N-1} \\ w_{N}\end{array}\right], \quad\) and \(\quad \mathbf{b}=\left[\begin{array}{c}-h^{2} r\left(x_{1}\right)+\left(1+\frac{h}{2} p\left(x_{1}\right)\right) w_{0} \\ -h^{2} r\left(x_{2}\right) \\ \vdots \\ -h^{2} r\left(x_{N-1}\right) \\ -h^{2} r\left(x_{N}\right)+\left(1-\frac{h}{2} p\left(x_{N}\right)\right) w_{N+1}\end{array}\right]\).

SPL tersebut akan diselesaikan dengan metode faktorisasi Crout (lihat algoritma 6.7). (basicly ini nyari inverse A secara linear, makanya runtime dari algortima ini adalah \(O(n)\))

Algoritma dari metode beda hingga linear:

function [xt,w]=linfdm(p,q,r,a_boundary,b_boundary,alpha,beta,n)
  h=(b_boundary-a_boundary)/(n+1); %stepsize
  a=zeros(n,1); %diagonal sistem persamaannya
  b=zeros(n,1); % right diagonal sistem persamaannya
  c=zeros(n,1); %left diagonal sistem persamaannya
  d=zeros(n,1); %vektor b (Ay=b) pada sistem persamaannya
  l=zeros(n,1); % main diagonal of lower triangle matrix
  u=zeros(n,1); %right diagonal of upper triangle matrix
  z= zeros(n,1); %solution of Lz=b
  w=zeros(n+1,1); %solusi aproksimasi dengan linear fdm
  xt=[a_boundary:h:b_boundary]; %mesh_point
  x=a_boundary+h;

  %konstruksi matrix tridiagonalnya
  a(1)=2+(h^2)*q(x);
  b(1)= -1+(h/2)*p(x);
  d(1)=-h^2*r(x) +(1+(h/2)*p(x))*alpha;

  for i = 2:n-1
    x= a_boundary+i*h;
    a(i)=2+h^2*q(x); %diagonal
    b(i)=-1+(h/2)*p(x);
    c(i)=-1-(h/2)*p(x);
    d(i)=-h^2*r(x);
  endfor

  x=b_boundary-h;
  a(n)=2+h^2*q(x);
  c(n)=-1-(h/2)*p(x);
  d(n)=-h^2*r(x)+(1-(h/2)*p(x))*beta;

  %matriks tridiagonalnya sudah didapatkan,
  %akan diselesaikan dengan LU Decomposition (crout factorization)

  l(1)= a(1);
  u(1)=b(1)/a(1);
  z(1)=d(1)/l(1);

  for i= 2:n-1
    l(i)=a(i)-c(i)*u(i-1);
    u(i)=b(i)/l(i);
    z(i)=(d(i)-c(i)*z(i-1))/l(i);

  endfor

  l(n)=a(n)-c(n)*u(n-1);
  z(n)=(d(n)-c(n)*z(n-1))/l(n);

  %konstruksi akhir w-nya
  w(n+1)=beta;
  w(n)=z(n);
  for i = n-1:-1:1
    w(i)=z(i)-u(i)*w(i+1);
  endfor

  w=[alpha;w];
  xt=transpose(xt);

endfunction

Akan kita uji dengan masalah syarat batas:

\[ \begin{aligned} y^{\prime \prime} & =-\frac{4}{x} y^{\prime}-\frac{2}{x^2} y+\frac{2 \ln x}{x^2}, \quad 1 \leq x \leq 2 \\ y(1) & =\frac{1}{2}, \quad y(2)=\ln 2 \end{aligned} \] Solusi eksak: \[ y(x)=\frac{4}{x}-\frac{2}{x^2}+\ln x-\frac{3}{2} \]

p= @(x) (-4/x); %function p(x)
q= @(x) (-2/x^2);%function q(x)
r=@(x) 2*log(x)/(x^2); %function r(x)
a_boundary=1; %batas kiri domain
b_boundary=2; %batas kanan domain
n=19; %banyaknya partisi (agar h=0.05 pilih n=19)
alpha=0.5; %y(a)=alpha
beta=log(2); %y(b)=beta
[x_grid,w]=linfdm(p,q,r,a_boundary,b_boundary,alpha,beta,n) %memangil fungsinya

f_anal=@(x)4./x -2./(x.^2) +log(x)-1.5;
sol_anal=f_anal(x_grid)
error=abs(sol_anal-w);

[x_grid,w,sol_anal,error]


%bikint tabel dan grafiknya :D

fplot(f_anal, [a_boundary,b_boundary],'b')
hold on;
scatter(x_grid,w,'r')
legend('solusi analitik', 'solusi linear FDM');
legend("location", "northwest");

Metode ini digunakan untuk mengaproksimasi masalah linear dalam bentuk:

\[ \begin{gathered} y^{\prime \prime}=f\left(x, y, y^{\prime}\right), \quad a \leq x \leq b \\ y(a)=\alpha, y(b)=\beta \end{gathered} \]

Aproksimasi menggunakan metode ini serupa dengan saat menggunakan metode beda hingga linear, dengan perbedaan kita juga menambahkan metode Newton dalam penyelesaiannya.

Algoritma dari metode beda hingga nonlinear:

function [t_grid,w]=nonlinear_FDM_naive(f,f_y,f_yprime,a,b,n,alpha,beta,max_iter,TOL)
  h=(b-a)/(n+1); %sepsize
  w=zeros(n,1); %vektor solusi aproksimasi
  t_grid=[a:h:b]; %mesh_poitnya
  J=zeros(n,n); %matriks jacobian
  F=zeros(n,1); %vektor fungsi  F=(f_1,f_2,...,f_n) yang dievaluasi di x_k

  for i=1:n %inisialisasi solusi awal
    w(i)=alpha+i*(beta-alpha)/(b-a)*h;
  endfor
  k=1;
  while k<=max_iter %lakukan iterasi jika masih belum didapat kriteria stopnya

    %solve nonlinear sistem tersebut dengan metode newton
    x=a+h;
    %kontruksi matriks Jacobian, dan vektor F-nya
    t=(w(2)-alpha)/(2*h);
    J(1,1)=2+h^2*f_y(x,w(1),t); %main diagoanal
    J(1,2)=-1+(h/2)*f_yprime(x,w(1),t); %right diagonal
    F(1)=(2*w(1)-w(2)-alpha+h^2*f(x,w(1),t));
    for i =2:n-1
      x=a+i*h;
      t=(w(i+1)-w(i-1))/(2*h);
      J(i,i)=2+h^2*f_y(x,w(i),t); %main diagoanal
      J(i,i+1)=-1+(h/2)*f_yprime(x,w(i),t); %main diagoanal
      J(i,i-1)=-1-(h/2)*f_yprime(x,w(i),t); %left diagoanal
      F(i)=(2*w(i)-w(i+1)-w(i-1)+h^2*f(x,w(i),t));
    endfor
     x=b-h;
     t=(beta-w(n-1))/(2*h);
     J(n,n)=2+h^2*f_y(x,w(n),t); %main diagonal
     J(n,n-1)=-1-(h/2)*f_yprime(x,w(n),t); %right diagonal
     F(n)=(2*w(n)-w(n-1)-beta+h^2*f(x,w(n),t));



    v=inverse(J)*F; %vector v adalah product dari J^-1 F
    w= w-v; % lakukan update nilai pada w

    if norm(v,2)<= TOL %kriteria stop jika norm(v)<=toleransinya
      break;
     else
        k=k+1; %jika belum memenuhi kriteria stop terus lanjut iterasinya (memperbaiki nilai w)
    endif
  endwhile
  w=[alpha ; w ; beta]; %konstruksi akhir w
  t_grid=transpose(t_grid); % %transpose meshpoint
  % untuk konsistensi dimensi saja

endfunction

Gunakan metode beda hingga nonlinear dengan \(h=0.1\) dan toleransi \(10^{-4}\) untuk mengaproksimasi BVP berikut: \[ \begin{aligned} y^{\prime \prime} & =y^{\prime}+2(y-\ln x)^3-\frac{1}{x}, \quad 2 \leq x \leq 3 \\ y(2) & =\frac{1}{2}+\ln 2, \quad y(3)=\frac{1}{3}+\ln 3 \end{aligned} \] Solusi eksak: \[ y(x)=\frac{1}{x}+\ln x \]

f=@(x,y,yp) yp+2*(y-log(x))^3-1/x ; %fungsi f pada y=f(x,y,y')
f_y=@(x,y,yp) 6*(y-log(x))^2; %turunan fungsi f terhadap y
f_yp=@(x,y,yp) 1; %turunan fungsi f terhadap yprime
a=2; %left boundary
b=3; %right boundary
alpha=0.5+log(2); %y(a)
beta=1/3+ log(3); %y(b)
n=9; %banyaknya partisi (pilih n=9 sehingga h=0.1)
maxiter=30; %masksimal iterasi newton methodnya
TOL=10^(-4); %toleransi nilai (untuk kriteria stop)

%memanggil fungsi nonlinear_FDM_naive
[x_grid,w]=nonlinear_FDM_naive(f,f_y,f_yp,a,b,n,alpha,beta,maxiter,TOL)
f_anal= @(x) 1./x +log(x); %sol analitik

%membuat grafiknya
fplot(f_anal, [a,b],'b')
hold on;
scatter(x_grid,w,'r')
legend('solusi analitik', 'solusi linear FDM');
legend("location", "northwest");




%membuat tabel saja.

sol_anal=f_anal(x_grid); %sol analitik di meshpoint
error=abs(w-sol_anal); %error
[x_grid,w,sol_anal,error]