sciforge_lib/maths/pde/
wave.rs

1pub fn wave_equation_1d(u: &[f64], u_prev: &[f64], dx: f64, dt: f64, c: f64) -> Vec<f64> {
2    let n = u.len();
3    let r = c * c * dt * dt / (dx * dx);
4    let mut u_new = vec![0.0; n];
5    u_new[0] = u[0];
6    u_new[n - 1] = u[n - 1];
7    for i in 1..n - 1 {
8        u_new[i] = 2.0 * u[i] - u_prev[i] + r * (u[i + 1] - 2.0 * u[i] + u[i - 1]);
9    }
10    u_new
11}
12
13pub fn wave_initial_step(u0: &[f64], v0: &[f64], dx: f64, dt: f64, c: f64) -> Vec<f64> {
14    let n = u0.len();
15    let r = c * c * dt * dt / (dx * dx);
16    let mut u1 = vec![0.0; n];
17    u1[0] = u0[0];
18    u1[n - 1] = u0[n - 1];
19    for i in 1..n - 1 {
20        u1[i] = u0[i] + dt * v0[i] + 0.5 * r * (u0[i + 1] - 2.0 * u0[i] + u0[i - 1]);
21    }
22    u1
23}
24
25pub fn wave_equation_2d(
26    u: &[Vec<f64>],
27    u_prev: &[Vec<f64>],
28    dx: f64,
29    dy: f64,
30    dt: f64,
31    c: f64,
32) -> Vec<Vec<f64>> {
33    let ny = u.len();
34    let nx = u[0].len();
35    let rx = c * c * dt * dt / (dx * dx);
36    let ry = c * c * dt * dt / (dy * dy);
37    let mut u_new = vec![vec![0.0; nx]; ny];
38    for j in 1..ny - 1 {
39        for i in 1..nx - 1 {
40            u_new[j][i] = 2.0 * u[j][i] - u_prev[j][i]
41                + rx * (u[j][i + 1] - 2.0 * u[j][i] + u[j][i - 1])
42                + ry * (u[j + 1][i] - 2.0 * u[j][i] + u[j - 1][i]);
43        }
44    }
45    u_new
46}
47
48pub fn courant_number(c: f64, dt: f64, dx: f64) -> f64 {
49    c * dt / dx
50}
51
52pub fn dalembert_solution(x: f64, t: f64, c: f64, f_init: fn(f64) -> f64) -> f64 {
53    0.5 * (f_init(x - c * t) + f_init(x + c * t))
54}
55
56pub fn wave_energy_density(
57    u: &[f64],
58    u_prev: &[f64],
59    dx: f64,
60    dt: f64,
61    c: f64,
62    rho: f64,
63) -> Vec<f64> {
64    let n = u.len();
65    let mut energy = vec![0.0; n];
66    for i in 1..n - 1 {
67        let du_dt = (u[i] - u_prev[i]) / dt;
68        let du_dx = (u[i + 1] - u[i - 1]) / (2.0 * dx);
69        energy[i] = 0.5 * rho * (du_dt * du_dt + c * c * du_dx * du_dx);
70    }
71    energy
72}
73
74pub fn absorbing_boundary(u: &mut [f64], u_prev: &[f64], dx: f64, dt: f64, c: f64) {
75    let n = u.len();
76    let alpha = c * dt / dx;
77    u[0] = u_prev[1] + (alpha - 1.0) / (alpha + 1.0) * (u[1] - u_prev[0]);
78    u[n - 1] = u_prev[n - 2] + (alpha - 1.0) / (alpha + 1.0) * (u[n - 2] - u_prev[n - 1]);
79}
80
81pub fn string_vibration_mode(x: f64, length: f64, mode: u32) -> f64 {
82    (mode as f64 * std::f64::consts::PI * x / length).sin()
83}
84
85pub fn wave_equation_1d_implicit(u: &[f64], u_prev: &[f64], dx: f64, dt: f64, c: f64) -> Vec<f64> {
86    let n = u.len();
87    let r = c * c * dt * dt / (dx * dx);
88    let mut a_diag = vec![0.0; n];
89    let mut b_diag = vec![0.0; n];
90    let mut c_diag = vec![0.0; n];
91    let mut d = vec![0.0; n];
92    b_diag[0] = 1.0;
93    d[0] = u[0];
94    b_diag[n - 1] = 1.0;
95    d[n - 1] = u[n - 1];
96    for i in 1..n - 1 {
97        a_diag[i] = -r;
98        b_diag[i] = 1.0 + 2.0 * r;
99        c_diag[i] = -r;
100        d[i] = 2.0 * u[i] - u_prev[i];
101    }
102    wave_thomas(&a_diag, &b_diag, &c_diag, &mut d);
103    d
104}
105
106pub fn cfl_check(c: f64, dt: f64, dx: f64) -> bool {
107    c * dt / dx <= 1.0
108}
109
110pub fn wave_reflection_coefficient(z1: f64, z2: f64) -> f64 {
111    (z2 - z1) / (z2 + z1)
112}
113
114pub fn wave_transmission_coefficient(z1: f64, z2: f64) -> f64 {
115    2.0 * z2 / (z2 + z1)
116}
117
118pub fn standing_wave(x: f64, t: f64, k: f64, omega: f64, amplitude: f64) -> f64 {
119    2.0 * amplitude * (k * x).sin() * (omega * t).cos()
120}
121
122pub fn wave_phase_velocity(omega: f64, k: f64) -> f64 {
123    omega / k
124}
125
126pub fn wave_group_velocity(domega: f64, dk: f64) -> f64 {
127    domega / dk
128}
129
130pub fn wave_total_energy(u: &[f64], u_prev: &[f64], dx: f64, dt: f64, c: f64, rho: f64) -> f64 {
131    let n = u.len();
132    let mut total = 0.0;
133    for i in 1..n - 1 {
134        let du_dt = (u[i] - u_prev[i]) / dt;
135        let du_dx = (u[i + 1] - u[i - 1]) / (2.0 * dx);
136        total += 0.5 * rho * (du_dt * du_dt + c * c * du_dx * du_dx) * dx;
137    }
138    total
139}
140
141pub fn spherical_wave_amplitude(r: f64, amplitude_0: f64, r0: f64) -> f64 {
142    if r < 1e-30 {
143        return 0.0;
144    }
145    amplitude_0 * r0 / r
146}
147
148pub fn wave_packet_gaussian(x: f64, t: f64, k0: f64, sigma: f64, c: f64) -> f64 {
149    let xi = x - c * t;
150    (-xi * xi / (4.0 * sigma * sigma)).exp() * (k0 * xi).cos()
151}
152
153pub fn wave_superposition(
154    x: f64,
155    t: f64,
156    amplitudes: &[f64],
157    frequencies: &[f64],
158    wave_numbers: &[f64],
159) -> f64 {
160    let mut val = 0.0;
161    for i in 0..amplitudes
162        .len()
163        .min(frequencies.len())
164        .min(wave_numbers.len())
165    {
166        val += amplitudes[i] * (wave_numbers[i] * x - frequencies[i] * t).sin();
167    }
168    val
169}
170
171pub fn wave_impedance(density: f64, velocity: f64) -> f64 {
172    density * velocity
173}
174
175fn wave_thomas(a: &[f64], b: &[f64], c: &[f64], d: &mut [f64]) {
176    let n = d.len();
177    let mut cp = vec![0.0; n];
178    let mut dp = d.to_vec();
179    cp[0] = c[0] / b[0];
180    dp[0] = d[0] / b[0];
181    for i in 1..n {
182        let m = b[i] - a[i] * cp[i - 1];
183        cp[i] = if i < n - 1 { c[i] / m } else { 0.0 };
184        dp[i] = (d[i] - a[i] * dp[i - 1]) / m;
185    }
186    d[n - 1] = dp[n - 1];
187    for i in (0..n - 1).rev() {
188        d[i] = dp[i] - cp[i] * d[i + 1];
189    }
190}
sciforge_lib/maths/pde/wave.rs

sciforge_lib/maths/pde/
wave.rs
