# PRAKTIKUM 13
`Solusi Persamaan Differensial Biasa (PDB) 2`

1. Runge-Kutta Fehlberg (RKF45)
2. Runge-Kutta untuk Sistem Persamaan Differensial Biasa
3. Persamaan Differensial Biasa (PDB) ordo tinggi (lebih dari 1)
4. PDB dengan Masalah Nilai Batas
    1. Metode _Linear-Shooting_
    2. Metode Beda-Hingga (_Finite-Difference_)


<hr style="border:2px solid black"> </hr>

# 1 Runge-Kutta Fehlberg (RKF45)
RKF45 merupakan modifikasi metode RK4 untuk memperoleh metode dengan ukuran langkah h yang adaptif, artinya ukuran langkah akan menyesuaikan bentuk fungsi (h kecil jika fungsi curam, h besar jika fungsi landai).
Solusi pada setiap iterasinya dihitung sebanyak dua kali, masing-masing menggunakan metode berorde 4 dan metode berorde 5.
Selisih kedua solusi pada setiap iterasinya digunakan untuk menentukan ukuran langkah h pada iterasi berikutnya. Ukuran langkah h dapat membesar atau mengecil sesuai kebutuhan. 


```julia
#%%METODE RUNGE-KUTTA FEHLBERG 4/5
#%
#% Digunakan untuk mencari solusi persamaan differensial
#%   dy/dt = f(t,y) dengan masalah nilai awal y(a) = y0
#%
#% sol = rkf45(f,a,b,y0,M,delta)
#% Input  : f    -> fungsi f(t,y)
#%          a,b  -> batas bawah dan atas solusi MNA
#%          y0   -> nilai awal y(a)=y0
#%          M    -> banyaknya sub-interval awal
#%          delta-> toleransi per langkah yang diberikan
#% Output : sol  -> solusi PD, sol=[T,Y]
#%
#% Digunakan Sebagai Pedoman Praktikum Metode Numerik
#%
#% Lihat juga : heun, taylor, rungekutta 
function rkf45(f,a,b,y0,M,delta)
  M = round(M)
  a2=1/4;b2=1/4;a3=3/8;b3=3/32;c3=9/32;
  a4=12/13;b4=1932/2197;c4=-7200/2197;d4=7296/2197;
  a5=1;b5=439/216;c5=-8;d5=3680/513;e5=-845/4104;
  a6=1/2;b6=-8/27;c6=2;d6=-3544/2565;e6=1859/4104;f6=-11/40;
  r1=1/360;r3=-128/4275;r4=-2197/75240;r5=1/50;r6=2/55;
  n1=25/216;n3=1408/2565;n4=2197/4104;n5=-1/5;
  big=1e15;
  h=(b-a)/M;
  hmin=h/64; 
  hmax=h*64; 
  maxi=200; 
  j=1;
  Y = y0;
  T = a;
  br= b-0.001*abs(b);
  err=NaN
  while T[j]<b
    if (T[j]+h)>br;h=b-T[j];end
    #% Hitung koefisien
    k1=h*f(T[j],Y[j]);
    y2=Y[j]+b2*k1;
    k2=h*f(T[j]+a2*h,y2);
    y3=Y[j]+b3*k1+c3*k2;
    k3=h*f(T[j]+a3*h,y3);
    y4=Y[j]+b4*k1+c4*k2+d4*k3;
    k4=h*f(T[j]+a4*h,y4);
    y5=Y[j]+b5*k1+c5*k2+d5*k3+e5*k4;
    k5=h*f(T[j]+a5*h,y5);
    y6=Y[j]+b6*k1+c6*k2+d6*k3+e6*k4+f6*k5;
    k6=h*f(T[j]+a6*h,y6);
    err=abs(r1*k1+r3*k3+r4*k4+r5*k5+r6*k6);
    ynew=Y[j]+n1*k1+n3*k3+n4*k4+n5*k5;
    #% Perbarui ukuran langkah
    if (err<delta) || (h<2*hmin) 
      Y = [Y; ynew];
      if (T[j]+h)>br
        T = [T; b];
      else
        T = [T; T[j]+h];
      end
      j=j+1;
    end
    if (err==0)
      s=0;
    else
      s=0.84*(delta*h/err)^(0.25);
    end
    if (s<0.75)&&(h>2*hmin);h = h/2;end
    if (s>1.50)&&(2*h<hmax);h = 2*h;end
    if abs(Y[j])>big || maxi==j;break;end
  end
  sol = [T Y];
  return sol
end
```

### Contoh 1
Diberikan persamaan differensial biasa $$y'=1+y^2$$ dengan masalah nilai awal $y(0)=0$. Berikut merupakan solusi numerik dari PDB tersebut menggunakan metode RKF45 pada $ t\in [0,1.4]$ dengan toleransi kesalahan $2\times10^{-5}$ dan ukuran langkah awal $h_0=0.2$.


```julia
# Hitung Solusi
f(t,y) = 1+y^2
a = 0
b = 1.4
y0 = 0
h0 = 0.2
M  = (b-a)/h0
delta = 2*10^-5
sol = rkf45(f,a,b,y0,M,delta)
```


```julia
using Plots
```


```julia
# Plot Solusi
t = sol[:,1];
y = sol[:,2];
plot(t,y,legend = :false)
scatter!(t,y)
```

# 2 Runge-Kutta untuk Sistem PDB
Sistem persamaan differensial merupakan masalah yang sering ditemui dalam pemodelan matematika. Sistem persamaan differensial biasanya digambarkan dalam suatu bidang fase dan bidang solusi yang dapat diselesaikan secara numerik menggunakan metode Runge-Kutta orde-4 seperti berikut.


```julia
#%%METODE RUNGE-KUTTA ORDE 4 untuk sistem PD
#%
#% Digunakan untuk mencari solusi sistem persamaan 
#%  differensial dx/dt = f(t,x,y) dengan x(a) = x0
#%               dy/dt = g(t,x,y) dengan y(a) = y0
#%
#% sol = rungekuttasistem(f,a,b,y0,M)
#% Input  : f    -> fungsi berisi f(t,[x,y]) dan g(t,[x,y])
#%          a,b  -> batas bawah dan atas solusi MNA
#%          y0   -> nilai awal y(a)=[x0,y0]
#%          M    -> banyaknya sub-interval
#% Output : sol  -> solusi PD, sol=[T,X,Y]
#%
#% Digunakan Sebagai Pedoman Praktikum Metode Numerik
#%
#% Lihat juga : heun, taylor, rungekutta
function rungekuttasistem(f,a,b,y0,M)
  M = Int(M)
  h = (b-a)/M;
  T = a:h:b;
  Y = Array{Float64}(undef,M+1,length(y0))
  Y[1,:] = y0;
  for k = 1:M
    f1 = f(T[k]     ,Y[k,:]        );
    f2 = f(T[k]+h/2 ,Y[k,:]+f1*h/2 );
    f3 = f(T[k]+h/2 ,Y[k,:]+f2*h/2 );
    f4 = f(T[k]+h   ,Y[k,:]+f3*h   );
    Y[k+1,:] = Y[k,:] + h/6*(f1+2*f2+2*f3+f4);
  end
  sol = [T Y];
end
```

### Contoh 2
Diberikan sistem persamaan differensial seperti berikut. 
\begin{align}  
\dot{x}= x+2y \\
\dot{y}=3x+2y 
\end{align}
dengan $ x(0)=6 $ dan $ y(0)=4 $.


Berikut merupakan langkah-langkah untuk menggambarkan plot bidang solusi dari sistem persamaan differensial di atas pada interval $ t\in[0, 0.2] $ dengan ukuran langkah $ h = 0.02 $.


```julia
# Hitung Solusi, misalkan x = z(1) dan y = z(2)
f(t,z) = [ z[1]+2*z[2] , 3*z[1]+2*z[2] ]
a = 0
b = 0.2
y0 = [6, 4]
h = 0.02
M = (b-a)/h
sol = rungekuttasistem(f,a,b,y0,M)
```


```julia
# Plot Fase
t = sol[:,1];
x = sol[:,2];
y = sol[:,3];
plt = plot(x,y,label="Bidang Fase",legend=:topleft,xlabel="x",ylabel="y")
```


```julia
# Plot Solusi
t = sol[:,1];
x = sol[:,2];
y = sol[:,3];
plt = plot(t,x,label="solusi x",legend=:topleft)
plot!(t,y,label="solusi y")
```

# 3. Persamaan Differensial Biasa (PDB) ordo tinggi (lebih dari 1)
Selain persamaan differensial orde pertama, sering kali dijumpai masalah matematika dalam bentuk persamaan differensial orde ke-2 atau lebih tinggi. Masalah persamaan differensial orde ke-2 atau lebih tinggi dapat diselesaikan dengan cara mentransformasi persamaan tersebut ke dalam bentuk sistem persamaan differensial orde pertama, kemudian diselesaikan menggunakan metode Runge-Kutta orde-4.

### Contoh 3
Diketahui suatu persamaan differensial orde-2 seperti berikut.
\begin{equation}\label{eq:13 orde2}
x''(t)+4x'(t)+5x(t)=0
\end{equation}
dengan nilai awal $ x(0)=3 $ dan $ x'(0)=-5 $. Berikut merupakan penyelesaian masalah tersebut menggunakan metode Runge-Kutta pada interval $ [0,5] $ dengan $ 50 $ sub-interval dan perbandingan solusi numerik terhadap solusi eksaknya, yaitu $$  x(t)=3e^{-2t}\cos(t)+e^{-2t}\sin(t)  $$

#### Transformasi persamaan differensial orde-2 menjadi sistem persamaan differensial orde-1

Misalkan, $ x'(t)=y(t) $ sehingga persamaan differensial memiliki bentuk 
\begin{align*}
				&\  x''(t)+4x'(t)+5x(t)=0 \\
\Leftrightarrow &\ y'(t)+4y(t)+5x(t)=0 \\
\Leftrightarrow &\ y'(t)=-4y(t)-5x(t)
\end{align*}
dan nilai awal $ x'(0)=-5 $ menjadi $ y(0)=-5 $. 

Secara lengkap, sistem baru yang equivalen dengan persamaan differensial di atas dapat dituliskan sebagai berikut.
\begin{align} 
x'&=y \\
y'&=-4y-5x						
\end{align}
dengan $x(0)=3$ dan $ y(0)=-5 $ 


```julia
# Hitung Solusi
f(t,z) = [ z[2] , -4*z[2]-5*z[1] ]
a = 0
b = 5
y0 = [3, -5]
M  = 50
sol = rungekuttasistem(f,a,b,y0,M)
```


```julia
# Plot Solusi
t = sol[:,1];
x = sol[:,2];
plt = plot(t,x,label="x(t)")
```


```julia
# Hitung Galat
xt(t) = 3*exp(-2*t)*cos(t)+exp(-2*t)*sin(t)
galat = abs.(x .- xt.(t))
```


```julia
# Plot Galat
plt = plot(t,galat)
```

# 4 PDB dengan Masalah Nilai Batas
Masalah nilai batas merupakan bentuk lain dari masalah persamaan differensial biasa. 

Masalah nilai batas memiliki bentuk umum berupa persamaan differensial orde 2 yaitu 

$$ x''=f(x',x,t) $$ 

dengan nilai batas pada selang $ [a,b] $, yaitu $ x(a)=\alpha $ dan $ x(b)=\beta $. 

Pada materi ini, akan dipelajari dua metode untuk mencari solusi numerik masalah nilai batas, yaitu metode \textit{linear shooting} dan \textit{finite-difference}.

## A. Metode _Linear-Shooting_
Ide dasar dari metode _linear shooting_ adalah melakukan transformasi masalah nilai batas menjadi dua persamaan differensial orde-2 dengan masalah nilai awal.

Misalkan, diketahui masalah nilai batas $ x''= p(t)x'+q(t)x+r(t) $ dengan nilai batas $ x(a)=\alpha $ dan $ x(b)=\beta $. Transformasi dari masalah nilai batas tersebut adalah

$ u''=p(t)u'+q(t)u+r(t) $ dengan $ u(a)=\alpha $ dan $ u'(a)=0 $

$ v''=p(t)v'+q(t)v      $ dengan $ v(a)=0 $ dan $ v'(a)=1 $

Persamaan tersebut dapat diselesaikan menggunakan metode Runge-Kutta orde-4 dengan mengubahnya menjadi bentuk sistem persamaan differensial, sehingga akan didapatkan solusi $ u(t) $ dan $ v(t) $. Selanjutnya, solusi masalah nilai batas $ x(t) $ dapat dicari dengan Persamaan berikut.

\begin{equation}\label{eq:13 linshoot}
x(t) = u(t)+\dfrac{\beta-u(b)}{v(b)}v(t)
\end{equation}


```julia
#%%METODE LINEAR SHOOTING
#%
#% Digunakan untuk mencari solusi persamaan differensial
#%       x''=p(t)x'+q(t)x+r(t) 
#% dengan masalah nilai batas x(a) = alpha , x(b)=beta
#%
#% fungsi ini membutuhkan adanya fungsi rungekuttasistem
#%
#% solusi = linearshooting(F1,F2,a,b,alpha,beta,M)
#% Input  : F1,F2     -> SPD hasil transformasi u dan v
#%          a,b       -> batas bawah dan atas solusi MNA
#%          alpha,beta-> nilai awal x(a)=alpha,x(b)=beta
#%          M         -> banyaknya sub-interval
#% Output : solusi    -> solusi PD, sol=[T,X]
#%
#% Digunakan Sebagai Pedoman Praktikum Metode Numerik
#%
#% Lihat juga : rkf45, findiff, rungekuttasistem 
function linearshooting(F1,F2,a,b,alpha,beta,M)
  M = Int(M)
  Za = [alpha,0];
  sol = rungekuttasistem(F1,a,b,Za,M);
  U = sol[:,2];
  
  Za = [0,1];
  sol = rungekuttasistem(F2,a,b,Za,M);
  V = sol[:,2];
  
  T = sol[:,1];
  X = U + (beta-U[M+1])*V/V[M+1];
  solusi = [T X];  
end
```

### Contoh 4
Diketahui persamaan masalah nilai batas seperti berikut.

\begin{equation}\label{eq:13 kasus 13.3}
x''(t)=\dfrac{2t}{1+t^2}x'(t)-\dfrac{2}{1+t^2}x(t)+1  
\end{equation}

dengan nilai batas $ x(0)=1.25 $ dan $ x(4)=-0.95 $ pada interval $ [0,4] $.

Berikut merupakan langkah-langkah untuk menyelesaikan masalah nilai batas tersebut menggunakan metode _linear shooting_.

#### Transformasi masalah nilai batas menjadi masalah nilai awal orde-2.

Berdasarkan Persamaan transformasi linear shooting, dua masalah nilai awal orde-2 yang setara dengan masalah nilai batas di atas adalah

$u''=\dfrac{2t}{1+t^2}u'-\dfrac{2}{1+t^2}u+1$ dengan $ u(0)=1.25 $ dan $ u'(0)=0 $

Serta,

$v''=\dfrac{2t}{1+t^2}v'-\dfrac{2}{1+t^2}v$ dengan $ v(0)=0 $ dan $ v'(0)=1 $.


#### Transformasi masing-masing masalah nilai awal orde-2 menjadi sistem persamaan differensial orde-1.
Masing-masing masalah nilai awal orde-2 di atas dapat ditulis dalam bentuk sistem persamaan linear orde pertama, yaitu 

$u'  =u_2 $

$u_2'=\dfrac{2t}{1+t^2}u_2-\dfrac{2}{1+t^2}u+1$

dengan $ u(0)=1.25 $ dan $ u_2(0)=0 $. Serta,

$v'	 =v_2$

$v_2' =\dfrac{2t}{1+t^2}v_2-\dfrac{2}{1+t^2}v$

dengan $ v(0)=0 $ dan $ v_2(0)=1 $.


```julia
# Hitung Solusi
F1(t,z) = [ z[2], 2*t/(1+t^2)*z[2]-2/(1+t^2)*z[1]+1 ]
F2(t,z) = [ z[2], 2*t/(1+t^2)*z[2]-2/(1+t^2)*z[1] ]
a = 0
b = 4
alpha = 1.25
beta = -0.95
M = 40
solusi = linearshooting(F1,F2,a,b,alpha,beta,M)
```


```julia
# Plot Solusi
t = solusi[:,1];
x = solusi[:,2];
plt = plot(t,x,legend=:false)
```

# B. Metode Beda-Hingga (Finite-Difference)
Ide dasar dari metode _finite-difference_ adalah mengubah persamaan differensial dengan masalah nilai batas menjadi formula beda-hingga. Bentuk SPL dari beda-hingga dapat dilihat pada bahan perkuliahan


```julia
#%%METODE FINITE-DIFFERENCE
#
#% Digunakan untuk mencari solusi persamaan differensial
#%       x''=p(t)x'+q(t)x+r(t) 
#% dengan masalah nilai batas x(a) = alpha , x(b)=beta
#
#% solusi = findiff(p,q,r,a,b,alpha,beta,M)
# Input  : p,q,r     -> fungsi p(t) q(t) dan r(t)
#%          a,b       -> batas bawah dan atas solusi MNA
#%          alpha,beta-> nilai awal x(a)=alpha,x(b)=beta
#          M         -> banyaknya sub-interval
#% Output : solusi    -> solusi PD, sol=[T,X]
#%
#% Digunakan Sebagai Pedoman Praktikum Metode Numerik
#%
#% Lihat juga : rkf45, linearshooting, rungekuttasistem 
using LinearAlgebra
function findiff(p,q,r,a,b,alpha,beta,M)
  h = (b-a)/M;
  T = a:h:b;
  T = T[2:end-1];
  #% Bangun Matriks B
  B = -h^2*r.(T);
  B[1] = B[1] + (1+h/2*p(T[1]))*alpha;
  B[end] = B[end] + (1-h/2*p(T[end]))*beta;
  # Bangun Matriks A - Bagian Diagonal 
  Ad = 2 .+h^2*q.(T);
  #% Bangun Matriks A - Bagian Bawah Diagonal 
  Tbawah = T[2:end];
  Abawah = -1 .-h/2*p.(Tbawah);
  #% Bangun Matriks A - Bagian Atas Diagonal 
  Tatas = T[1:end-1];
  Aatas = -1 .+h/2*p.(Tatas);
  A = Tridiagonal(Abawah,Ad,Aatas)
  #% Selesaikan AX=B
  X = A\B;
  T = [a; T; b];
  X = [alpha; X ;beta];
  solusi = [T X];
end

```

### Contoh 5
Diketahui persamaan masalah nilai batas seperti berikut.

\begin{equation}
x''(t)=\dfrac{2t}{1+t^2}x'(t)-\dfrac{2}{1+t^2}x(t)+1  
\end{equation}

dengan nilai batas $ x(0)=1.25 $ dan $ x(4)=-0.95 $ pada interval $ [0,4] $.

Berikut merupakan langkah-langkah untuk menyelesaikan masalah nilai batas tersebut menggunakan metode _finite-difference_.


```julia
# Hitung Solusi
p(t)= 2*t/(1+t^2);
q(t)= -2/(1+t^2);
r(t)= 1 + 0*t;
a = 0
b = 4
alpha = 1.25
beta = -0.95
M = 40
solusi = findiff(p,q,r,a,b,alpha,beta,M)
```


```julia
# Plot Solusi
t = solusi[:,1];
x = solusi[:,2];
plt = plot(t,x,legend=:false)
```

<hr style="border:2px solid black"> </hr>

# Soal Latihan
Kerjakan soal berikut pada saat kegiatan praktikum berlangsung.

`Nama: ________`

`NIM: ________`

### Soal 1
Diberikan persamaan differensial

$$ \dfrac{dy}{dt} = -y-2t-1 ,\ \text{ dengan } y(0)=2 $$

Gunakan metode Runge-Kutta-Fehlberg 4/5  (RKF45) untuk menyelesaikan PD di atas dengan nilai toleransi $ \delta=10^{-6} $, $ 10^{-12} $ dan $ 10^{-16} $, kemudian bandingkan dengan solusi eksaknya yaitu $ y(t)=e^{-t}-2t+1 $ dan gambarkan solusi beserta titik-titik pembagian ukuran langkah seperti pada **Contoh 1**.


```julia

```

### Soal 2
Ulangi langkah-langkah pada **Contoh 2** untuk menggambarkan plot bidang fase dan bidang solusi dari sistem persamaan berikut.

$$
\begin{align} 
	\dot{x}&= -2x-y \\ 
	\dot{y}&=   x-y \\ 
\end{align}
$$

dengan $ x(0)=6 $ dan $ y(0)=0 $ pada interval $ t\in[0,\ 5] $ dengan ukuran langkah $ h=0.1 $.


```julia

```

### Soal 3
Diketahui suatu persamaan differensial orde-2 seperti berikut.

$$ 2x''(t)-5x'(t)-3x(t)=45e^{2t} $$

dengan nilai awal $ x(0)=2 $ dan $ x'(0)=1 $. Selesaikan masalah tersebut menggunakan Runge-Kutta pada interval $ [0,\ 2] $ dengan $ h=0.05 $, kemudian bandingkan solusi numerik terhadap solusi eksaknya yaitu $ x(t)=4e^{-t/2}+7e^{3t}-9e^{2t} $ dengan langkah-langkah seperti pada **Contoh 3**.


```julia

```

### Soal 4
Diketahui persamaan masalah nilai batas seperti berikut.

$$ x''(t)=-\frac{2}{t}\ x'(t)+\frac{2}{t^2}\ x(t)+\frac{\sin(t)}{t^2} $$

dengan nilai batas $ x(1)=-0.02 $ dan $ x(6)=0.02 $ pada interval $ [1,6] $. Gunakan metode _linear shooting_ untuk menyelesaikan masalah nilai batas di atas dengan langkah-langkah pada **Contoh 4**.


```julia

```

### Soal 5
Diketahui persamaan masalah nilai batas seperti berikut.

$$ x''(t)=-\frac{2}{t}\ x'(t)+\frac{2}{t^2}\ x(t)+\frac{\sin(t)}{t^2} $$

dengan nilai batas $ x(1)=-0.02 $ dan $ x(6)=0.02 $ pada interval $ [1,6] $. Gunakan metode _finite difference_ untuk menyelesaikan masalah nilai batas di atas dengan langkah-langkah pada **Contoh 5**.


```julia

```

54b4b418ab811dce36b9064e611b855bcc8fb410
ipynb
Jupyter Notebook
notebookpraktikum/Praktikum 13.ipynb
mkhoirun-najiboi/metnum.jl
a6e35d04dc277318e32256f9b432264157e9b8f4
MIT
notebookpraktikum/Praktikum 13.ipynb
mkhoirun-najiboi/metnum.jl
a6e35d04dc277318e32256f9b432264157e9b8f4
MIT
notebookpraktikum/Praktikum 13.ipynb
mkhoirun-najiboi/metnum.jl
a6e35d04dc277318e32256f9b432264157e9b8f4
MIT
Qwen/Qwen-72B
1. YES
2. YES