Persamaan diubah menjadi persamaan diferensial $$\frac{∂^2 P}{dt^2 }=c^2 \frac{∂^2 P_0}{dx^2 }-2α \frac{∂P}{dt} $$

Diskritisasi kedua ruang dan waktu memberikan:

$$\frac{P_{i,j+1}-2P_{i,j}+P_{i,j-1}}{\Delta t^2}=c^2\frac{P_{i+1,j}-2P_{i,j}+P_{i-1,j}}{\Delta x^2}-2\alpha\frac{P_{i,j+1}-P_{i,j-1}}{2\Delta t}$$

Dengan menggunakan notasi indeks Crank-Nicolson, kita dapat menyusun persamaan sebagai berikut:

$$[P_{i,j+1}=2rP_{i+1,j}+\left(1-2r\right)P_{i,j}+rP_{i-1,j}-\frac{2\alpha\Delta t}{2+\alpha\Delta t}\left(P_{i,j+1}-P_{i,j-1}\right)]$$

di mana $$(r=\frac{c^2\Delta t^2}{\Delta x^2})$$

 1import numpy as np
 2import matplotlib.pyplot as plt
 3import math
 4
 5def sound_wave_amplitude(x, P0, alpha, rho, c, f, t):
 6    omega = 2 * np.pi * f
 7    k = omega / c
 8    lambd = 2 * np.pi / k
 9    
10    amplitude = P0 * np.exp(-alpha * x) * np.sin((2 * np.pi / lambd) * x - omega * t)
11    
12    return amplitude
13
14def simulate_damped_wave_crank_nicolson(P0, alpha, rho, c, f, L, T, Nx, Nt):
15    dx = L / (Nx - 1)
16    dt = T / (Nt - 1)
17    
18    x_values = np.linspace(0, L, Nx)
19    t_values = np.linspace(0, T, Nt)
20    
21    u = np.zeros((Nx, Nt))
22    
23    u[:, 0] = sound_wave_amplitude(x_values, P0, alpha, rho, c, f, 0)
24    
25    r = alpha * dt / (2 * dx)
26    for n in range(0, Nt - 1):
27        A = np.eye(Nx) + r * np.diag(np.ones(Nx - 1), 1) - r * np.diag(np.ones(Nx - 1), -1)
28        B = np.eye(Nx) - r * np.diag(np.ones(Nx - 1), 1) + r * np.diag(np.ones(Nx - 1), -1)
29        
30        u[:, n + 1] = np.linalg.solve(A, np.dot(B, u[:, n]))
31    
32    return x_values, t_values, u
33
34# Parameters
35P0 = 1.0
36alpha = 0.01
37rho = 1000.0
38L = 100.0
39T = 20
40Nx = 100  # Increase spatial grid points
41Nt = 100  # Increase time grid points
42
43# Bulk modulus and frequency
44B = 1.96 * 10**9
45C = (B * 1) / rho
46c = math.sqrt(C)
47f = 100.0
48
49# Solve with Crank-Nicolson method
50x_values, t_values, u = simulate_damped_wave_crank_nicolson(P0, alpha, rho, c, f, L, T, Nx, Nt)
51
52# Plot the result
53plt.plot(x_values, sound_wave_amplitude(x_values, P0, alpha, rho, c, f, T), label='Solusi Analitik', linestyle='--')
54plt.plot(x_values, u[:, -1], label='Solusi numerik')
55plt.xlabel('Jarak (m)')
56plt.ylabel('Tekanan Suara (Pa)')
57plt.title('Perbandingan solusi numerik dan Analitik')
58plt.legend()
59plt.savefig('perbandingan solusi analitik dan numerik.jpg', dpi=300)
60plt.show()
61
62# Calculate error at a specific time (e.g., at the final time)
63t_final_index = 1
64error = np.abs(u[:, t_final_index] - sound_wave_amplitude(x_values, P0, alpha, rho, c, f, T))
65
66print('Average error: {}'.format(np.mean(error, axis=0)))
67x_values_target = 100
68index_x_target = np.argmin(np.abs(x_values - x_values_target))
69print('Amplitude at x={} at t={}: {}'.format(x_values_target, t_final_index, u[index_x_target, t_final_index]))
70print('Analytical Amplitude at x={} at t={}: {}'.format(x_values_target, T, sound_wave_amplitude(x_values_target, P0, alpha, rho, c, f, T)))
71print("c: {}".format(c))
72print("k: {}".format(2 * np.pi * f / c))