Modul 6 Persamaan Diferensial Numerik: Metode Shooting dan Linear Finite Difference untuk Masalah Nilai Batas PDB

Modul 6 Persamaan Diferensial Numerik: Metode Shooting dan Linear Finite Difference untuk Masalah Nilai Batas PDB

Kembali ke Persamaan Diferensial Numerik

Di modul ini, kita akan membahas beberapa metode untuk masalah nilai batas untuk PDB, yaitu:

1. Linear Shooting
2. Nonlinear Shotting
3. Linear Finite Difference

Review: Runge-Kutta orde 4 untuk sistem

Shooting method untuk masalah nilai batas melibatkan sistem persamaan diferensial. Kita akan menggunakan kode metode Runge-Kutta orde 4 untuk sistem dari modul sebelumnya:

Function file rko4_sysm.m - nama file harus sama dengan nama fungsi

function [t, w] = rko4_sysm(cell_f, a, b, N, alphas)
  m = length(cell_f);

  h = (b - a) / N;
  t = zeros(N + 1, 1);
  w = zeros(m, N + 1);
  t(1) = a;
  w(:, 1) = alphas;

  k1 = zeros(m, 1);
  k2 = zeros(m, 1);
  k3 = zeros(m, 1);
  k4 = zeros(m, 1);
  for i = 1 : N
    t(i + 1) = t(i) + h;

    for j = 1 : m
      k1(j) = h * cell_f{j}(t(i), w(:, i));
    endfor

    for j = 1 : m
      k2(j) = h * cell_f{j}(t(i) + (h / 2), w(:, i) + (k1 / 2));
    endfor

    for j = 1 : m
      k3(j) = h * cell_f{j}(t(i) + (h / 2), w(:, i) + (k2 / 2));
    endfor

    for j = 1 : m
      k4(j) = h * cell_f{j}(t(i + 1), w(:, i) + k3);
    endfor

    for j = 1 : m
      w(j, i + 1) = w(j, i) + (k1(j) + 2 * k2(j) + 2 * k3(j) + k4(j)) / 6;
    endfor
  endfor
endfunction

Sebenarnya tidak harus metode Runge-Kutta orde 4 untuk sistem. Boleh ditukar dengan metode lainnya untuk sistem, misalnya metode Adams predictor-corrector orde 4 untuk sistem.

Linear Shooting

Bentuk umum, ide utama, penyederhanaan

Linear Shooting merupakan metode untuk menyelesaikan sejenis masalah nilai batas untuk PDB, yaitu yang berbentuk:

\(y'' = f\left(x,y,y'\right) = p(x)y' + q(x)y + r(x), \;a\leq x\leq b\)

\(y(a)=\alpha, \;y(b)=\beta\)

dengan

• \(p(x), q(x), r(x)\) adalah fungsi kontinu dalam \(x\)
• \(q(x) > 0\) pada \([a,b]\) agar dijamin ada solusi unik

Cara penyelesaiannya:

1. Selesaikan MNA PDB orde 2 berikut, solusinya disebut \(y_1 \left(x\right)\):

\[y'' = p(x)y' + q(x)y + r(x), \quad a \le x \le b, \quad y\left(a\right) = \alpha, \quad y'\left(a\right) = 0\]

2. Selesaikan MNA PDB orde 2 berikut, solusinya disebut \(y_2 \left(x\right)\)

\[y'' = p(x)y' + q(x)y, \quad a \le x \le b, \quad y\left(a\right) = 0, \quad y'\left(a\right) = 1\]

3. Solusi akhirnya adalah

\[y\left(x\right) = y_1 \left(x\right) + \frac{\beta - y_1 \left(b\right)}{y_2 \left(b\right)} y_2 \left(x\right)\]

Kita bisa menuliskan kedua MNA PDB orde 2 tersebut masing-masing sebagai sistem PDB orde 1, seperti biasa dengan permisalan \(u_1(x) = y(x)\) dan \(u_2(x) = y'(x)\).

Sehingga, langkahnya menjadi:

1. Selesaikan sistem PDB orde 1 berikut. Kemudian solusi \(u_1(x)\) disebut \(y_1(x)\) dan solusi \(u_2(x)\) disebut \(y_1'(x)\).

\[\begin{aligned} u_1'(x) &= u_2(x) \\ u_2'(x) &= p(x) u_2(x) + q(x) u_1(x) + r(x) \\ u_1 (a) &= \alpha, \quad u_2 (a) = 0 \end{aligned}\]

2. Selesaikan sistem PDB orde 1 berikut. Kemudian solusi \(u_1(x)\) disebut \(y_2(x)\) dan solusi \(u_2(x)\) disebut \(y_2'(x)\).

\[\begin{aligned} u_1'(x) &= u_2(x) \\ u_2'(x) &= p(x)u_2(x) + q(x)u_1(x) \\ u_1 (a) &= 0, \quad u_2(a) = 1 \end{aligned}\]

3. Solusi akhirnya adalah

\[y(x) = y_1(x) + \frac{\beta - y_1(b)}{y_2(b)} y_2(x)\]

Kalau perlu,

\[y'(x) = y_1'(x) + \frac{\beta - y_1(b)}{y_2(b)} y_2'(x)\]

Function file (dari pseudocode di buku)

Function file linshoot_pseudocode.m - nama file harus sama dengan nama fungsi

function [x_i, w_1i, w_2i] = linshoot_pseudocode(p, q, r, a, b, n, alpha, beta)
  h = (b - a)/n;
  u = [alpha ; 0];
  v = [0 ; 1];
  x_i = w_1i = w_2i = [];
  for i = 1:n
    x = a + (i-1)*h;

    k_11 = h * u(2,i);
    k_12 = h * (p(x)*u(2,i) + q(x)*u(1,i) + r(x));

    k_21 = h * (u(2,i)+(k_12/2));
    k_22 = h * (p(x+(h/2))*(u(2,i)+(k_12/2)) + q(x+(h/2))*(u(1,i)+(k_11/2)) + r(x+(h/2)));

    k_31 = h * (u(2,i)+(k_22/2));
    k_32 = h * (p(x+(h/2))*(u(2,i)+(k_22/2)) + q(x+(h/2))*(u(1,i)+(k_21/2)) + r(x+(h/2)));

    k_41 = h * (u(2,i)+k_32);
    k_42 = h * (p(x+h)*(u(2,i)+k_32) + q(x+h)*(u(1,i)+k_31) + r(x+h));

    u(1,i+1) = u(1,i) + ((k_11 + 2*k_21 + 2*k_31 + k_41)/6);
    u(2,i+1) = u(2,i) + ((k_12 + 2*k_22 + 2*k_32 + k_42)/6);

    kp_11 = h * v(2,i);
    kp_12 = h * (p(x)*v(2,i) + q(x)*v(1,i));

    kp_21 = h * (v(2,i) + (kp_12/2));
    kp_22 = h * (p(x+(h/2))*(v(2,i)+(kp_12/2)) + q(x+(h/2))*(v(1,i)+(kp_11/2)));

    kp_31 = h * (v(2,i)+(kp_22/2));
    kp_32 = h * (p(x+(h/2))*(v(2,i)+(kp_22/2)) + q(x+(h/2))*(v(1,i)+(kp_21/2)));

    kp_41 = h * (v(2,i)+kp_32);
    kp_42 = h * (p(x+h)*(v(2,i)+kp_32) + q(x+h)*(v(1,i)+kp_31));

    v(1,i+1) = v(1,i) + (kp_11 + 2*kp_21 + 2*kp_31 + kp_41)/6;
    v(2,i+1) = v(2,i) + (kp_12 + 2*kp_22 + 2*kp_32 + kp_42)/6;
  endfor

  w = [alpha ; ((beta - u(1,(n+1))) / v(1,(n+1)))];
  x_i(1) = a;
  w_1i(1) = w(1,1);
  w_2i(1) = w(2,1);

  for i = 2:(n+1)
    W1 = u(1,i) + w(2,1)*v(1,i);
    W2 = u(2,i) + w(2,1)*v(2,i);
    x = a + (i-1)*h;
    x_i(i) = x;
    w_1i(i) = W1;
    w_2i(i) = W2;
  endfor
endfunction

Function file (lebih sederhana)

Function file linear_shooting.m - nama file harus sama dengan nama fungsi

function [x, w1, w2] = linear_shooting(p, q, r, a, b, N, alph, bet)
  % sistem PDB yang pertama
  u1_aksen = @(x, u) u(2);
  u2_aksen = @(x, u) p(x)*u(2) + q(x)*u(1) + r(x);
  [x, w_pers1] = rko4_sysm({u1_aksen, u2_aksen}, a, b, N, [alph, 0]);
  y1_b = w_pers1(1, N+1);
  
  % sistem PDB yang kedua
  u1_aksen = @(x, u) u(2);
  u2_aksen = @(x, u) p(x)*u(2) + q(x)*u(1);
  [x, w_pers2] = rko4_sysm({u1_aksen, u2_aksen}, a, b, N, [0, 1]);
  y2_b = w_pers2(1, N+1);
  
  % solusi akhir (superposisi)
  w_akhir = w_pers1 + (bet - y1_b)/(y2_b) * w_pers2;
  % dipisah jadi w1,i (aproksimasi y(x)) dan w2,i (aproksimasi y'(x))
  w1 = w_akhir(1, :)'; % ditranspos agar menjadi vektor kolom
  w2 = w_akhir(2, :)';
endfunction

Contoh Linear Shooting

\(y'' = -\frac{2}{x}y' + \frac{2}{x^2}y + \frac{\sin(\ln(x))}{x^2}, \; 1\leq x\leq 2\)

\(y(1)=1,\; y(2)=2\)

dengan \(N=10\)

dan solusi eksak:

\(y(x)=c_1x+\frac{c_2}{x^2} - \frac{3}{10}\sin(\ln(x))-\frac{1}{10}\cos(\ln(x))\)

\(y'(x)=c_1 - \frac{2c_2}{x^3} - \frac{3\cos(\ln(x))}{10x} + \frac{\sin(\ln(x))}{10x}\)

\(c_2 = \frac{1}{70}(8-12\sin(\ln(2)) - 4\cos(\ln(2)))\)

\(c_1 = \frac{11}{10}-c_2\)

Berikut code script file untuk permasalahan di atas menggunakan metode linear shooting:

Script file coba_linear_shooting.m - nama file bebas

p = @(x) (-2 ./ x);
q = @(x) (2 ./ (x .^ 2));
r = @(x) (sin(log(x)) ./ (x .^ 2));
a = 1;
b = 2;
N = 10;
alph = 1;
bet = 2;

[xi, w1i, w2i] = linear_shooting(p, q, r, a, b, N, alph, bet);

% solusi eksak y(x)
c2 = (8-12*sin(log(2)) - 4*cos(log(2)))/70;
c1 = (11/10) - c2;
sln = @(x) (c1*x + (c2 ./ x.^2) - (3/10)*sin(log(x)) - (1/10)*cos(log(x)));
y_eksak = sln(xi);

% menghitung error
err_w1i = abs(w1i - y_eksak);
err1_total = sum(err_w1i);

% tampilkan
format long;
disp("Tabel aproksimasi w1,i, solusi eksak y(x), dan error:");
[xi, w1i, y_eksak, err_w1i]
disp("Error total (norm L1):");
disp(err1_total);
format;

% plot untuk y(x)
figure;
hold on;
fplot(sln, [a,b], 'k');
scatter(xi, w1i, 'r');
title("Aproksimasi y(x)");
legend("Eksak", "Aproksimasi (w1,i)");
legend('location', 'northwest');

Tabel aproksimasi w1,i, solusi eksak y(x), dan error:
ans =

   1.000000000000000   1.000000000000000   1.000000000000000                   0
   1.100000000000000   1.092629164133552   1.092629298481288   0.000000134347735
   1.200000000000000   1.187084706810955   1.187084840483685   0.000000133672730
   1.300000000000000   1.283382266283346   1.283382364079130   0.000000097795784
   1.400000000000000   1.381445891533503   1.381445951696987   0.000000060163484
   1.500000000000000   1.481159386366171   1.481159416999814   0.000000030633643
   1.600000000000001   1.582392449986370   1.582392460756381   0.000000010770011
   1.700000000000001   1.685013962277612   1.685013961734097   0.000000000543514
   1.800000000000001   1.788898539692082   1.788898534641947   0.000000005050134
   1.900000000000001   1.893929513621671   1.893929509211183   0.000000004410488
   2.000000000000001   2.000000000000000   2.000000000000000   0.000000000000000

Error total (norm L1):
4.773875239560965e-07

Jika kita run script file tersebut, maka program akan mengeluarkan dua macam output, yaitu tabel serta plot perbandingan solusi eksak dan aproksimasi dengan errornya.

Apabila diperlukan bukan hanya \(y(x)\) (yaitu \(w_{1,i}\)) tetapi juga \(y'(x)\) (yaitu \(w_{2,i}\)), kodenya menjadi seperti berikut:

Script file coba2_linear_shooting.m - nama file bebas

p = @(x) (-2 ./ x);
q = @(x) (2 ./ (x .^ 2));
r = @(x) (sin(log(x)) ./ (x .^ 2));
a = 1;
b = 2;
N = 10;
alph = 1;
bet = 2;

[xi, w1i, w2i] = linear_shooting(p, q, r, a, b, N, alph, bet);

% solusi eksak y(x) dan y'(x)
c2 = (8-12*sin(log(2)) - 4*cos(log(2)))/70;
c1 = (11/10) - c2;
sln = @(x) (c1*x + (c2 ./ x.^2) - (3/10)*sin(log(x)) - (1/10)*cos(log(x)));
sln_p = @(x) (c1 - (2*c2 ./ x.^3) - 3*cos(log(x))./(10*x) + sin(log(x))./(10*x));
y_eksak = sln(xi);
yp_eksak = sln_p(xi);

% menghitung error
err_w1i = abs(w1i - y_eksak);
err_w2i = abs(w2i - yp_eksak);
err1_total = sum(err_w1i);
err2_total = sum(err_w2i);

% tampilkan
format long;
disp("Tabel aproksimasi w1,i, solusi eksak y(x), dan error:");
[xi, w1i, y_eksak, err_w1i]
disp("Tabel aproksimasi w2,i, solusi eksak y'(x), dan error:");
[xi, w2i, yp_eksak, err_w2i]
disp("Error total (norm L1) untuk w1,i:");
disp(err1_total);
disp("Error total (norm L1) untuk w2,i:");
disp(err2_total);
format;

% plot untuk y(x)
figure;
hold on;
fplot(sln, [a,b], 'k');
scatter(xi, w1i, 'r');
title("Aproksimasi y(x)");
legend("Eksak", "Aproksimasi (w1,i)");
legend('location', 'northwest');

% plot untuk y'(x)
figure;
hold on;
fplot(sln_p, [a,b], 'k');
scatter(xi, w2i, 'b');
title("Aproksimasi y'(x)");
legend("Eksak", "Aproksimasi (w2,i)");
legend('location', 'northwest');

Tabel aproksimasi w1,i, solusi eksak y(x), dan error:
ans =

   1.000000000000000   1.000000000000000   1.000000000000000                   0
   1.100000000000000   1.092629164133552   1.092629298481288   0.000000134347735
   1.200000000000000   1.187084706810955   1.187084840483685   0.000000133672730
   1.300000000000000   1.283382266283346   1.283382364079130   0.000000097795784
   1.400000000000000   1.381445891533503   1.381445951696987   0.000000060163484
   1.500000000000000   1.481159386366171   1.481159416999814   0.000000030633643
   1.600000000000001   1.582392449986370   1.582392460756381   0.000000010770011
   1.700000000000001   1.685013962277612   1.685013961734097   0.000000000543514
   1.800000000000001   1.788898539692082   1.788898534641947   0.000000005050134
   1.900000000000001   1.893929513621671   1.893929509211183   0.000000004410488
   2.000000000000001   2.000000000000000   2.000000000000000   0.000000000000000

Tabel aproksimasi w2,i, solusi eksak y'(x), dan error:
ans =

 Columns 1 through 3:

   1.000000000000000e+00   9.176213963846825e-01   9.176210394808360e-01
   1.100000000000000e+00   9.352828620960821e-01   9.352826025486669e-01
   1.200000000000000e+00   9.538386694942848e-01   9.538385751339796e-01
   1.300000000000000e+00   9.719773228811912e-01   9.719773623999391e-01
   1.400000000000000e+00   9.890965253486431e-01   9.890966553369642e-01
   1.500000000000000e+00   1.004953219518007e+00   1.004953404763131e+00
   1.600000000000001e+00   1.019487696052374e+00   1.019487911769632e+00
   1.700000000000001e+00   1.032732443129011e+00   1.032732672951544e+00
   1.800000000000001e+00   1.044763943100940e+00   1.044764176648869e+00
   1.900000000000001e+00   1.055676941937622e+00   1.055677172867598e+00
   2.000000000000001e+00   1.065570770445955e+00   1.065570995061434e+00

 Column 4:

   3.569038465878194e-07
   2.595474152267130e-07
   9.436030512510740e-08
   3.951874794072552e-08
   1.299883211069996e-07
   1.852451241290964e-07
   2.157172589445366e-07
   2.298225325603198e-07
   2.335479285520137e-07
   2.309299755864913e-07
   2.246154791052390e-07

Error total (norm L1) untuk w1,i:
4.773875239560965e-07
Error total (norm L1) untuk w2,i:
2.200196934865062e-06

Nonlinear Shooting

Bentuk umum, ide utama

Nonlinear Shooting digunakan untuk menyelesaikan masalah nilai batas berbentuk:

\(y'' = f(x, y, y'), \; a\leq x \leq b\)

\(y(a)=\alpha, \; y(b)=\beta\)

dengan \(f\) boleh berupa fungsi linier maupun nonlinier

Cara penyelesaiannya:

1. Tentukan toleransi \(\varepsilon\), dan pilih tebakan awal \(t_0\) (yaitu \(t_k\) sebelum iterasi pertama, yaitu dengan \(k=0\)). Kalau bingung, disarankan

\[t_0 = \frac{\beta - \alpha}{b-a}\]

2. Selesaikan MNA PDB orde 2 berikut, misalkan solusinya disebut \(w(x,t_k)\):

\[y'' = f(x,y,y'), \quad a \le x \le b, \quad y(a) = \alpha, \quad y'(a) = t_k\]

3. Periksa apakah \(\left|w(b,t_k) - \beta\right| \le \varepsilon\).

   • Kalau iya, selesai; solusi akhirnya adalah \(y(x) = w(x,t_k)\).
   • Kalau tidak, peroleh tebakan baru untuk \(t_i\) (misalnya dengan metode secant atau metode Newton), lalu kembali ke langkah 2.

Seperti biasa, kita bisa misalkan \(u_1(x) = y(x)\) dan \(u_2(x) = y'(x)\) agar MNA PDB orde 2 menjadi sistem PDB orde 1.

Cara penyelesaiannya menjadi:

1. Tentukan toleransi \(\varepsilon\), dan pilih tebakan awal \(t_0\) (yaitu \(t_k\) sebelum iterasi pertama, yaitu dengan \(k=0\)). Kalau bingung, disarankan

\[t_0 = \frac{\beta - \alpha}{b-a}\]

2. Selesaikan sistem PDB orde 1 berikut. Kemudian \(u_1(x)\) disebut \(w(x,t_k)\) dan \(u_2(x)\) disebut \(w'(x,t_k)\).

\[\begin{aligned} u_1'(x) &= u_2(x) \\ u_2'(x) &= f(x,u_1,u_2) \\ y(a) &= \alpha, \quad y'(a) = t_k \end{aligned}\]

3. Periksa apakah \(\left|w(b,t_k) - \beta\right| \le \varepsilon\).

   • Kalau iya, selesai; solusi akhirnya adalah \(y(x) = w(x,t_k)\).
   • Kalau tidak, peroleh tebakan baru untuk \(t_k\) (misalnya dengan metode secant atau metode Newton), lalu kembali ke langkah 2.

Function file (metode secant)

Function file nonlinear_shooting_secant.m - nama file harus sama dengan nama fungsi

function [x, w1, w2] = nonlinear_shooting_secant(f, a, b, N, alph, bet, tol, t0, t1)
  u1_aksen = @(x, u) u(2);
  u2_aksen = @(x, u) f(x, u(1), u(2));
  
  t_k_min_2 = t0;
  t_k_min_1 = t1;
  [x, w_k_min_2] = rko4_sysm({u1_aksen, u2_aksen}, a, b, N, [alph, t_k_min_2]);
  [x, w_k_min_1] = rko4_sysm({u1_aksen, u2_aksen}, a, b, N, [alph, t_k_min_1]);
  w_k = w_k_min_1;
  err = abs(w_k(1, N+1) - bet);
  while !(err <= tol)
    pembilang = (w_k_min_1(1,N+1) - bet) * (t_k_min_1 - t_k_min_2);
    penyebut = w_k_min_1(1,N+1) - w_k_min_2(1,N+1);
    t_k = t_k_min_1 - pembilang/penyebut;

    t_k_min_2 = t_k_min_1
    t_k_min_1 = t_k

    [x, w_k] = rko4_sysm({u1_aksen, u2_aksen}, a, b, N, [alph, t_k]);
    err = abs(w_k(1, N+1) - bet);
  endwhile
  % keluar loop artinya toleransi sudah terpenuhi
  
  % memisahkan w_k menjadi w1i dan w2i
  w1 = w_k(1, :)'; % transpos juga agar menjadi vektor kolom
  w2 = w_k(2, :)';
endfunction

Modifikasi untuk metode Newton

Untuk menggunakan metode Newton, diperlukan tidak hanya \(y(b,t)\) tetapi juga turunannya \(\frac{\partial y(b,t)}{\partial t}\) yang sayangnya tidak dimiliki.

Setelah penjabaran yang panjang di buku, ternyata bisa dimisalkan

\[z(x,t) = \frac{\partial y(x,t)}{\partial t}\]

dan nilai fungsi \(z\) ini ternyata bisa diperoleh dengan menyelesaikan suatu MNA PDB orde 2 (lagi). Sehingga, di tiap iterasi, ada dua MNA PDB orde 2 yang harus diselesaikan.

Langkah nonlinear shooting dengan metode Newton bisa ditulis:

1. Hitung rumus \(\frac{\partial f}{\partial y}(x,y,y')\) dan rumus \(\frac{\partial f}{\partial y'}(x,y,y')\) secara analitik.

2. Tentukan toleransi \(\varepsilon\), dan pilih tebakan awal \(t_0\) (yaitu \(t_k\) sebelum iterasi pertama, yaitu dengan \(k=0\)). Kalau bingung, disarankan

\[t_0 = \frac{\beta - \alpha}{b-a}\]

3. Selesaikan MNA PDB orde 2 berikut, misalkan solusinya disebut \(w(x,t_k)\):

\[y'' = f(x,y,y'), \quad a \le x \le b, \quad y(a) = \alpha, \quad y'(a) = t_k\]

4. Selesaikan MNA PDB orde 2 berikut:

\[z'' = \frac{\partial f}{\partial y}(x,y,y')z(x) + \frac{\partial f}{\partial y'}(x,y,y')z'(x), \quad z(a) = 0, \quad z'(a) = 1\]

5. Periksa apakah \(\left|w(b,t_k) - \beta\right| \le \varepsilon\).

   • Kalau iya, selesai; solusi akhirnya adalah \(y(x) = w(x,t_k)\).
   • Kalau tidak, kembali ke langkah 3 setelah memperoleh tebakan baru untuk \(t_k\): \[t_k = t_{k-1} - \frac{w(b, t_{k-1}) - \beta}{z(b, t_{k-1})}\]

Dengan permisalan \(u_1\) dan \(u_2\) agar PDB orde 2 menjadi sistem PDB orde 1, langkah-langkahnya menjadi:

1. Hitung rumus \(\frac{\partial f}{\partial y}(x,y,y')\) dan rumus \(\frac{\partial f}{\partial y'}(x,y,y')\) secara analitik.

2. Tentukan toleransi \(\varepsilon\), dan pilih tebakan awal \(t_0\) (yaitu \(t_k\) sebelum iterasi pertama, yaitu dengan \(k=0\)). Kalau bingung, disarankan

\[t_0 = \frac{\beta - \alpha}{b-a}\]

3. Selesaikan sistem PDB orde 1 berikut. Kemudian \(u_1(x)\) disebut \(w(x,t_k)\) dan \(u_2(x)\) disebut \(w'(x,t_k)\).

\[\begin{aligned} u_1'(x) &= u_2(x) \\ u_2'(x) &= f(x,u_1,u_2) \\ y(a) &= \alpha, \quad y'(a) = t_k \end{aligned}\]

4. Selesaikan sistem PDB orde 1 berikut. Kemudian \(u_1(x)\) disebut \(z(x,t_k)\) dan \(u_2(x)\) disebut \(z'(x,t_k)\).

   \[\begin{aligned} u_1'(x) &= u_2(x) \\ u_2'(x) &= \frac{\partial f}{\partial y}(x,y,y')u_1(x) + \frac{\partial f}{\partial y'}(x,y,y')u_2(x) \\ u_1(a) &= 0, \quad u_2(a) = 1 \end{aligned}\]

   Note: nilai \(y\) dan \(y'\) sebenarnya tergantung \(x\), sehingga sebaiknya ditulis \(y(x)\) dan \(y'(x)\):

   \[\begin{aligned} u_1'(x) &= u_2(x) \\ u_2'(x) &= \frac{\partial f}{\partial y}(x,y(x),y'(x))u_1(x) + \frac{\partial f}{\partial y'}(x,y(x),y'(x))u_2(x) \\ u_1(a) &= 0, \quad u_2(a) = 1 \end{aligned}\]

   Dalam perhitungan, nilai \(y(x)\) dan \(y'(x)\) bisa kita peroleh dari \(w(x,t_k)\) dan \(w'(x,t_k)\), yaitu nilai \(w_{1,i}\) dan \(w_{2,i}\) dari sistem yang sebelumnya.

5. Periksa apakah \(\left|w(b,t_k) - \beta\right| \le \varepsilon\).

   • Kalau iya, selesai; solusi akhirnya adalah \(y(x) = w(x,t_k)\).
   • Kalau tidak, kembali ke langkah 3 setelah memperoleh tebakan baru untuk \(t_k\): \[t_k = t_{k-1} - \frac{w(b, t_{k-1}) - \beta}{z(b, t_{k-1})}\]

Function file (dari pseudocode)

Function file nonlinshoot_pseudocode.m - nama file harus sama dengan nama fungsi

function [x_i, w_1i, w_2i] = nonlinshoot_pseudocode(f, fy, fyp, a, b, n, alpha, beta, m, tol)
  % m adalah maksimum iterasi

  h = (b - a)/n;
  k = 1;
  tk = (beta - alpha)/(b - a);
  x_i = w_1i = w_2i = [];
  while k <= m
    w = [alpha;tk];
    u = [0,1];
    for i = 1:n
      x = a + (i-1)*h;

      k_11 = h*w(2,i);
      k_12 = h*f(x, w(1,i), w(2,i));

      k_21 = h*(w(2,i)+(k_12/2));
      k_22 = h*f((x+(h/2)), (w(1,i)+(k_11/2)), (w(2,i)+(k_12/2)));

      k_31 = h*(w(2,i)+(k_22/2));
      k_32 = h*f((x+(h/2)), (w(1,i)+(k_21/2)), (w(2,i)+(k_22/2)));

      k_41 = h*(w(2,i)+k_32);
      k_42 = h*f((x+h), (w(1,i)+k_31), (w(2,i)+k_32));

      w(1,i+1) = w(1,i) + ((k_11 + 2*k_21 + 2*k_31 + k_41)/6);
      w(2,i+1) = w(2,i) + ((k_12 + 2*k_22 + 2*k_32 + k_42)/6);

      kp_11 = h*u(2);
      kp_12 = h*(fy(x, w(1,i), w(2,i))*u(1) + fyp(x, w(1,i), w(2,i))*u(2));

      kp_21 = h*(u(2) + (kp_12/2));
      kp_22 = h*(fy((x+(h/2)), w(1,i), w(2,i))*u(1) + fyp((x+(h/2)), w(1,i), w(2,i))*(u(2) + (kp_12/2)));

      kp_31 = h*(u(2)+(kp_22/2));
      kp_32 = h*(fy((x+(h/2)), w(1,i), w(2,i))*(u(1) + (kp_21/2)) + fyp((x+(h/2)), w(1,i), w(2,i))*(u(2) + (kp_22/2)));

      kp_41 = h*(u(2)+kp_32);
      kp_42 = h*(fy((x+h), w(1,i), w(2,i))*(u(1)+kp_31) + fyp((x+h), w(1,i), w(2,i))*(u(2) + kp_32));

      u(1) = u(1) + (kp_11 + 2*kp_21 + 2*kp_31 + kp_41)/6;
      u(2) = u(2) + (kp_12 + 2*kp_22 + 2*kp_32 + kp_42)/6;
    endfor

  if abs(w(1,n+1) - beta) <= tol       % jika sudah mencapai batas toleransi maka program berhenti
    for i = 1:(n+1)
      x = a+(i-1)*h;
      x_i(i) = x;
      w_1i(i) = w(1,i);
      w_2i(i) = w(2,i);
    endfor
    return
  endif
  tk = tk-((w(1,n+1) - beta)/u(1));
  k = k + 1;
  endwhile
  disp('max iteration')
endfunction

Function file (metode Newton)

Function file nonlinear_shooting_newton.m - nama file harus sama dengan nama fungsi

function [x_arr, w1i, w2i] = nonlinear_shooting_newton(f, fy, fyp, a, b, N, alph, bet, tol, M, t0)
  % kalau input t0 bukan angka, dianggap tidak memilih tebakan awal
  if isnumeric(t0)
    t_k = t0;
  else
    t_k = (bet-alph)/(b-a);
  endif

  % banyaknya iterasi
  k = 1;

  err = tol + 1;
  % selama belum memenuhi toleransi ataupun mencapai batas iterasi
  while (!(err <= tol) && k != M+1)
    % selesaikan sistem pertama
    u1_aksen = @(x, u) u(2);
    u2_aksen = @(x, u) f(x, u(1), u(2));
    [x_arr, w_sys] = rko4_sysm({u1_aksen, u2_aksen}, a, b, N, [alph, t_k]);

    % y(x) adalah w1,i dengan i adalah indeks dari x (perlu dicari)
    % carinya bisa dengan memilih nilai terakhir (indeks end) di x_arr yang <= x
    % (seandainya ada misalnya t_i + h/2, nilai w yang digunakan tetap di t_i)
    y = @(x) w_sys(1, find(x_arr <= x)(end));

    % y'(x) adalah w2,i dengan i adalah indeks dari x (perlu dicari)
    yp = @(x) w_sys(2, find(x_arr <= x)(end));

    % selesaikan sistem kedua
    u1_aksen = @(x, u) u(2);
    u2_aksen = @(x, u) fy(x, y(x), yp(x))*u(1) + fyp(x, y(x), yp(x))*u(2);
    [x_arr, z_sys] = rko4_sysm({u1_aksen, u2_aksen}, a, b, N, [0, 1]);

    % periksa toleransi, update t_k
    err = abs(w_sys(1, N+1) - bet);
    if !(err <= tol)
      t_k = t_k - (w_sys(1, N+1) - bet)/(z_sys(1, N+1));
    endif

    % lanjut iterasi selanjutnya
    k += 1;
  endwhile
  % keluar loop artinya toleransi sudah terpenuhi atau maks iterasi tercapai
  if (k == M)
    printf("Maks iterasi (%d) tercapai\n", M);
  endif

  % pisahkan w_sys menjadi w1i dan w2i
  w1i = w_sys(1, :)'; % transpos juga agar menjadi vektor kolom
  w2i = w_sys(2, :)';
endfunction

Contoh Nonlinear Shooting

\(y'' = \frac{1}{8}(32+2x^3-yy'), \; 1\leq x \leq 3\)

\(y(1) = 17, \; y(3)=43/3\)

dengan \(N=20\) dan toleransi \(=10^{-5}\)

dan solusi eksak:

\(y(x)=x^2 + \frac{16}{x}\)

\(y'(x)=2x - \frac{16}{x^2}\)

Hint:

\[\begin{aligned} f(x,y,y') &= \frac{1}{8}(32+2x^3-yy') \\ \frac{\partial f}{\partial y}(x,y,y') &= -\frac{1}{8}y' \\ \frac{\partial f}{\partial y'}(x,y,y') &= -\frac{1}{8}y \end{aligned}\]

Berikut code script file untuk permasalahan di atas menggunakan metode nonlinear shooting (dengan metode Newton):

Script file coba_nonlinear_shooting_newton.m - nama file bebas

f = @(x, y, yp) ((1/8)*(32 + 2 * x.^3 - y .* yp));
fy = @(x, y, yp) (-yp/8);
fyp = @(x, y, yp) (-y/8);
a = 1;
b = 3;
N = 20;
alph = 17;
bet = 43/3;
tol = 10^(-5);
M = -1; % maks iterasi. Nilai negatif artinya tidak ada maks

[xi, w1i, w2i] = nonlinear_shooting_newton(f, fy, fyp, a, b, N, alph, bet, tol, M, "");

% solusi eksak y(x)
sln = @(x) ((x .^ 2) + (16 ./ x));
y_eksak = sln(xi);

% menghitung error
err_w1i = abs(w1i - y_eksak);
err1_total = sum(err_w1i);

% tampilkan
format long;
disp("Tabel aproksimasi w1,i, solusi eksak y(x), dan error:");
[xi, w1i, y_eksak, err_w1i]
disp("Error total (norm L1):");
disp(err1_total);
format;

figure;
hold on;
fplot(sln, [a,b], 'k');
scatter(xi, w1i, 'r');
title("Aproksimasi y(x)");
legend("Eksak", "Aproksimasi (w1,i)");
legend('location', 'northeast');

Tabel aproksimasi w1,i, solusi eksak y(x), dan error:
ans =

 Columns 1 through 3:

   1.000000000000000e+00   1.700000000000000e+01   1.700000000000000e+01
   1.100000000000000e+00   1.575549578599484e+01   1.575545454545455e+01
   1.200000000000000e+00   1.477339050038047e+01   1.477333333333333e+01
   1.300000000000000e+00   1.399775337003275e+01   1.399769230769231e+01
   1.400000000000000e+00   1.338863063808904e+01   1.338857142857143e+01
   1.500000000000000e+00   1.291672136654978e+01   1.291666666666667e+01
   1.600000000000001e+00   1.256004909500447e+01   1.256000000000000e+01
   1.700000000000001e+00   1.230180789967695e+01   1.230176470588235e+01
   1.800000000000001e+00   1.212892629211319e+01   1.212888888888889e+01
   1.900000000000001e+00   1.203108455596395e+01   1.203105263157895e+01
   2.000000000000001e+00   1.200002684877419e+01   1.200000000000000e+01
   2.100000000000001e+00   1.202906982837228e+01   1.202904761904762e+01
   2.200000000000001e+00   1.211274528064585e+01   1.211272727272727e+01
   2.300000000000001e+00   1.224653596950516e+01   1.224652173913044e+01
   2.400000000000001e+00   1.242667752129748e+01   1.242666666666667e+01
   2.500000000000001e+00   1.265000785517754e+01   1.265000000000000e+01
   2.600000000000001e+00   1.291385135925930e+01   1.291384615384616e+01
   2.700000000000002e+00   1.321592880477185e+01   1.321592592592593e+01
   2.800000000000002e+00   1.355428656395277e+01   1.355428571428572e+01
   2.900000000000002e+00   1.392724047230069e+01   1.392724137931035e+01
   3.000000000000002e+00   1.433333091826060e+01   1.433333333333334e+01

 Column 4:

                       0
   4.124054029475133e-05
   5.716704713520926e-05
   6.106234043912195e-05
   5.920951761773097e-05
   5.469988311368468e-05
   4.909500447425330e-05
   4.319379459793993e-05
   3.740322430267895e-05
   3.192438500754236e-05
   2.684877418701603e-05
   2.220932465668568e-05
   1.800791857675677e-05
   1.423037472392252e-05
   1.085463081196281e-05
   7.855177535986968e-06
   5.205413140529913e-06
   2.878845918985462e-06
   8.496670478308488e-07
   9.070096620433787e-07
   2.415072742678603e-06

Error total (norm L1):
5.472579459873117e-04

Jika kita run script file tersebut, maka program akan mengeluarkan dua macam output, yaitu tabel serta plot perbandingan solusi eksak dan aproksimasi dengan errornya.

Apabila diperlukan bukan hanya \(y(x)\) (yaitu \(w_{1,i}\)) tetapi juga \(y'(x)\) (yaitu \(w_{2,i}\)), kodenya menjadi seperti berikut:

f = @(x, y, yp) ((1/8)*(32 + 2 * x.^3 - y .* yp));
fy = @(x, y, yp) (-yp/8);
fyp = @(x, y, yp) (-y/8);
a = 1;
b = 3;
N = 20;
alph = 17;
bet = 43/3;
tol = 10^(-5);
M = -1; % maks iterasi. Nilai negatif artinya tidak ada maks

[xi, w1i, w2i] = nonlinear_shooting_newton(f, fy, fyp, a, b, N, alph, bet, tol, M, "");

% solusi eksak y(x) dan y'(x)
sln = @(x) ((x .^ 2) + (16 ./ x));
sln_p = @(x) (2*x - 16 ./ (x.^2));
y_eksak = sln(xi);
yp_eksak = sln_p(xi);

% menghitung error
err_w1i = abs(w1i - y_eksak);
err_w2i = abs(w2i - yp_eksak);
err1_total = sum(err_w1i);
err2_total = sum(err_w2i);

% tampilkan
format long;
disp("Tabel aproksimasi w1,i, solusi eksak y(x), dan error:");
[xi, w1i, y_eksak, err_w1i]
disp("Tabel aproksimasi w2,i, solusi eksak y'(x), dan error:");
[xi, w2i, yp_eksak, err_w2i]
disp("Error total (norm L1) untuk w1,i:");
disp(err1_total);
disp("Error total (norm L1) untuk w2,i:");
disp(err2_total);
format;

figure;
hold on;
fplot(sln, [a,b], 'k');
scatter(xi, w1i, 'r');
title("Aproksimasi y(x)");
legend("Eksak", "Aproksimasi (w1,i)");
legend('location', 'northeast');

figure;
hold on;
fplot(sln_p, [a,b], 'k');
scatter(xi, w2i, 'b');
title("Aproksimasi y'(x)");
legend("Eksak", "Aproksimasi (w2,i)");
legend('location', 'northwest');

Tabel aproksimasi w1,i, solusi eksak y(x), dan error:
ans =

 Columns 1 through 3:

   1.000000000000000e+00   1.700000000000000e+01   1.700000000000000e+01
   1.100000000000000e+00   1.575549578599484e+01   1.575545454545455e+01
   1.200000000000000e+00   1.477339050038047e+01   1.477333333333333e+01
   1.300000000000000e+00   1.399775337003275e+01   1.399769230769231e+01
   1.400000000000000e+00   1.338863063808904e+01   1.338857142857143e+01
   1.500000000000000e+00   1.291672136654978e+01   1.291666666666667e+01
   1.600000000000001e+00   1.256004909500447e+01   1.256000000000000e+01
   1.700000000000001e+00   1.230180789967695e+01   1.230176470588235e+01
   1.800000000000001e+00   1.212892629211319e+01   1.212888888888889e+01
   1.900000000000001e+00   1.203108455596395e+01   1.203105263157895e+01
   2.000000000000001e+00   1.200002684877419e+01   1.200000000000000e+01
   2.100000000000001e+00   1.202906982837228e+01   1.202904761904762e+01
   2.200000000000001e+00   1.211274528064585e+01   1.211272727272727e+01
   2.300000000000001e+00   1.224653596950516e+01   1.224652173913044e+01
   2.400000000000001e+00   1.242667752129748e+01   1.242666666666667e+01
   2.500000000000001e+00   1.265000785517754e+01   1.265000000000000e+01
   2.600000000000001e+00   1.291385135925930e+01   1.291384615384616e+01
   2.700000000000002e+00   1.321592880477185e+01   1.321592592592593e+01
   2.800000000000002e+00   1.355428656395277e+01   1.355428571428572e+01
   2.900000000000002e+00   1.392724047230069e+01   1.392724137931035e+01
   3.000000000000002e+00   1.433333091826060e+01   1.433333333333334e+01

 Column 4:

                       0
   4.124054029475133e-05
   5.716704713520926e-05
   6.106234043912195e-05
   5.920951761773097e-05
   5.469988311368468e-05
   4.909500447425330e-05
   4.319379459793993e-05
   3.740322430267895e-05
   3.192438500754236e-05
   2.684877418701603e-05
   2.220932465668568e-05
   1.800791857675677e-05
   1.423037472392252e-05
   1.085463081196281e-05
   7.855177535986968e-06
   5.205413140529913e-06
   2.878845918985462e-06
   8.496670478308488e-07
   9.070096620433787e-07
   2.415072742678603e-06

Tabel aproksimasi w2,i, solusi eksak y'(x), dan error:
ans =

 Columns 1 through 3:

   1.000000000000000e+00  -1.400019602422796e+01  -1.400000000000000e+01
   1.100000000000000e+00  -1.102334186854095e+01  -1.102314049586777e+01
   1.200000000000000e+00  -8.711296684185600e+00  -8.711111111111109e+00
   1.300000000000000e+00  -6.867619909427612e+00  -6.867455621301771e+00
   1.400000000000000e+00  -5.363408505219970e+00  -5.363265306122443e+00
   1.500000000000000e+00  -4.111235293408942e+00  -4.111111111111106e+00
   1.600000000000001e+00  -3.050107650303584e+00  -3.049999999999994e+00
   1.700000000000001e+00  -2.136425654341337e+00  -2.136332179930791e+00
   1.800000000000001e+00  -1.338352959256321e+00  -1.338271604938266e+00
   1.900000000000001e+00  -6.322039283580119e-01  -6.321329639889148e-01
   2.000000000000001e+00  -6.200755926466517e-05   5.329070518200751e-15
   2.100000000000001e+00   5.718278547215568e-01   5.718820861678053e-01
   2.200000000000001e+00   1.094167447234676e+00   1.094214876033063e+00
   2.300000000000001e+00   1.575383898119827e+00   1.575425330812860e+00
   2.400000000000001e+00   2.022186112082093e+00   2.022222222222227e+00
   2.500000000000001e+00   2.439968644139854e+00   2.440000000000006e+00
   2.600000000000001e+00   2.833109007887408e+00   2.833136094674562e+00
   2.700000000000002e+00   3.205189382604409e+00   3.205212620027440e+00
   2.800000000000002e+00   3.559163917474358e+00   3.559183673469393e+00
   2.900000000000002e+00   3.897486371338108e+00   3.897502972651611e+00
   3.000000000000002e+00   4.222208482001776e+00   4.222222222222228e+00

 Column 4:

   1.960242279608337e-04
   2.013726731835419e-04
   1.855730744910744e-04
   1.642881258403506e-04
   1.431990975273578e-04
   1.241822978359508e-04
   1.076503035899457e-04
   9.347441054607941e-05
   8.135431805422755e-05
   7.096436909703741e-05
   6.200755926999424e-05
   5.423144624849829e-05
   4.742879838692815e-05
   4.143269303269470e-05
   3.611014013493730e-05
   3.135586015190484e-05
   2.708678715324098e-05
   2.323742303111942e-05
   1.975599503500902e-05
   1.660131350300631e-05
   1.374022045119716e-05

Error total (norm L1) untuk w1,i:
5.472579459873117e-04
Error total (norm L1) untuk w2,i:
1.741071134524930e-03

Linear Finite Difference

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}\]

Penyelesaiannya adalah dengan persamaan-persamaan berikut:

\[\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 disusun menjadi 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]\).

Menurut buku, SPL tersebut sebaiknya diselesaikan dengan metode faktorisasi Crout (algoritma 6.7). (Intinya, mumpung A adalah matriks tridiagonal, algoritma ini nyari inverse A secara linier, makanya runtime dari algortima ini adalah \(O(n)\))

Namun, tentunya kita bebas menyelesaikan SPLnya dengan cara apapun, misalnya dengan invers, OBE, atau bahkan dengan cara iteratif seperti Gauss-Seidel

Function file (dari pseudocode)

Function file linfdm_pseudocode.m - nama file harus sama dengan nama fungsi

function [xt,w]=linfdm_pseudocode(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

Function file (dengan solusi SPL secara langsung/invers)

Function file linear_fd_langsung.m - nama file harus sama dengan nama fungsi

function [x, w_grid] = linear_fd_langsung(p, q, r, a, b, N, alph, bet)
  % bikin array x
  h = (b - a) / (N+1);
  x = (a : h : b)'; % transpos juga agar menjadi vektor kolom

  % susun matriks A dan vektor b
  A = zeros(N, N);
  b = -h^2 * r(x(2:N+1));
  % kasus khusus untuk baris pertama
  b(1) += (1 + h/2 * p(x(2))) * alph;
  A(1, 1) += 2 + h^2 * q(x(2));
  A(1, 2) += -1 + h/2 * p(x(2));
  % kasus khusus untuk baris terakhir
  A(N, N-1) += -1 - h/2 * p(x(N+1));
  A(N, N) += 2 + h^2 * q(x(N+1));
  b(N) += (1 - h/2 * p(x(N+1))) * bet;
  % untuk baris kedua hingga kedua-terakhir
  for i = 2 : (N-1)
    A(i, i-1) += -1 - h/2 * p(x(i+1));
    A(i, i) += 2 + h^2 * q(x(i+1));
    A(i, i+1) += -1 + h/2 * p(x(i+1));
  endfor

  % selesaikan SPL
  w = A \ b;
  % w baru mengandung w1, ..., w_N

  % gabungkan dengan w0 (alpha) dan w_{N+1} (beta)
  w_grid = [alph w' bet]'; % transpos juga agar menjadi vektor kolom
endfunction

Contoh Linear Finite Difference

Akan kita uji dengan masalah nilai 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} \]

Script file coba_linear_fd_langsung.m - nama file bebas

p= @(x) (-4./x); %fungsi p(x)
q= @(x) (-2./x.^2); %fungsi q(x)
r=@(x) 2*log(x)./(x.^2); %fungsi r(x)
a=1; %batas kiri interval
b=2; %batas kanan interval
N=20; %banyaknya partisi
alph=0.5; %y(a)=alpha
bet=log(2); %y(b)=beta
[x_grid, w_grid] = linear_fd_langsung(p,q,r,a,b,N,alph,bet); %memangil fungsinya

% solusi eksak y(x) dan error
sln = @(x) 4./x - 2./(x.^2) + log(x) - 1.5;
y_eksak = sln(x_grid);
err = abs(y_eksak - w_grid);
err_total = sum(err); % norm L1 (taxicab/Manhattan)

% bikin tabel dan grafiknya :D

format long;
disp("Tabel aproksimasi, solusi eksak y(x), dan error:");
[x_grid, w_grid, y_eksak, err]
disp("Error total (norm L1):");
disp(err_total);
format;

hold on;
fplot(sln, [a,b], 'b');
scatter(x_grid, w_grid, 'r');
title("Aproksimasi y(x)");
legend("Eksak", "Aproksimasi");
legend('location', 'northwest');

Tabel aproksimasi, solusi eksak y(x), dan error:
ans =

   1.000000000000000   0.500000000000000   0.500000000000000                   0
   1.047619047619048   0.542451840408551   0.542387784229934   0.000064056178617
   1.095238095238095   0.575949844443050   0.575848904859791   0.000100939583258
   1.142857142857143   0.602401287601637   0.602281392624522   0.000119894977115
   1.190476190476190   0.623280489841670   0.623153387144777   0.000127102696893
   1.238095238095238   0.639736322556670   0.639609603256639   0.000126719300031
   1.285714285714286   0.652670546690794   0.652548996182141   0.000121550508653
   1.333333333333333   0.662795564639343   0.662682072451781   0.000113492187563
   1.380952380952381   0.670677452832402   0.670573630075180   0.000103822757222
   1.428571428571429   0.676768343081596   0.676674943938732   0.000093399142864
   1.476190476190476   0.681431010801427   0.681348221496374   0.000082789305052
   1.523809523809524   0.684957702680262   0.684885340076304   0.000072362603959
   1.571428571428571   0.687584665665992   0.687522313825702   0.000062351840290
   1.619047619047619   0.689503439735737   0.689450543640143   0.000052896095593
   1.666666666666667   0.690869694214655   0.690825623765991   0.000044070448665
   1.714285714285714   0.691810185166203   0.691774278510465   0.000035906655738
   1.761904761904762   0.692428265250871   0.692399857681941   0.000028407568931
   1.809523809523810   0.692808270884264   0.692786713692713   0.000021557191551
   1.857142857142857   0.693019033126842   0.693003705447643   0.000015327679199
   1.904761904761905   0.693116700585629   0.693107016390513   0.000009684195115
   1.952380952380952   0.693147019138584   0.693142430884514   0.000004588254070
   2.000000000000000   0.693147180559945   0.693147180559945   0.000000000000000

Error total (norm L1):
1.400919170378323e-03

Sayangnya, metode finite difference untuk masalah nilai batas tidak bisa menentukan aproksimasi \(y'(x)\).

Perhatikan errornya. Menurut buku Burden, untuk masalah nilai batas, metode shooting pada umumnya lebih akurat dibandingkan metode finite difference.