scirs2_integrate/pde/implicit/mod.rs
1//! Implicit methods for solving PDEs
2//!
3//! This module provides implementations of implicit time-stepping schemes for
4//! solving partial differential equations (PDEs). Implicit methods are generally
5//! more stable than explicit methods and can handle stiff problems more efficiently.
6//!
7//! ## Supported Methods
8//!
9//! * Crank-Nicolson method: A second-order implicit method that is A-stable
10//! * Backward Euler method: A first-order fully implicit method that is L-stable
11//! * Trapezoidal rule: A second-order implicit method
12//! * Alternating Direction Implicit (ADI) method: An efficient operator splitting method for
13//! multi-dimensional problems
14
15use scirs2_core::ndarray::{Array1, Array2, ArrayView1};
16use std::time::Instant;
17
18use crate::pde::finite_difference::FiniteDifferenceScheme;
19use crate::pde::{
20 BoundaryCondition, BoundaryConditionType, BoundaryLocation, Domain, PDEError, PDEResult,
21 PDESolution, PDESolverInfo,
22};
23
24/// Type alias for coefficient functions in PDEs
25type CoefficientFunction = Box<dyn Fn(f64, f64, f64) -> f64 + Send + Sync>;
26
27/// Type alias for initial condition function
28type InitialCondition = Box<dyn Fn(f64) -> f64 + Send + Sync>;
29
30// Re-export ADI implementation
31mod adi;
32pub use adi::{ADIResult, ADI2D};
33
34/// Available implicit time-stepping schemes
35#[derive(Debug, Clone, Copy, PartialEq)]
36pub enum ImplicitMethod {
37 /// Backward Euler method (first-order, L-stable)
38 BackwardEuler,
39
40 /// Crank-Nicolson method (second-order, A-stable)
41 CrankNicolson,
42
43 /// Trapezoidal rule (second-order, A-stable)
44 TrapezoidalRule,
45
46 /// Alternating Direction Implicit method (efficient for multi-dimensional problems)
47 ADI,
48}
49
50/// Options for implicit PDE solvers
51#[derive(Debug, Clone)]
52pub struct ImplicitOptions {
53 /// Implicit time-stepping method to use
54 pub method: ImplicitMethod,
55
56 /// Tolerance for iterative solvers
57 pub tolerance: f64,
58
59 /// Maximum number of iterations for linear system solver
60 pub max_iterations: usize,
61
62 /// Fixed time step (if None, adaptive time-stepping will be used)
63 pub dt: Option<f64>,
64
65 /// Minimum time step for adaptive time-stepping
66 pub min_dt: Option<f64>,
67
68 /// Maximum time step for adaptive time-stepping
69 pub max_dt: Option<f64>,
70
71 /// Time steps to save (if None, all steps will be saved)
72 pub save_every: Option<usize>,
73
74 /// Print detailed progress information
75 pub verbose: bool,
76}
77
78impl Default for ImplicitOptions {
79 fn default() -> Self {
80 ImplicitOptions {
81 method: ImplicitMethod::CrankNicolson,
82 tolerance: 1e-6,
83 max_iterations: 100,
84 dt: Some(0.01),
85 min_dt: Some(1e-6),
86 max_dt: Some(0.1),
87 save_every: None,
88 verbose: false,
89 }
90 }
91}
92
93/// Result of implicit method solution
94pub struct ImplicitResult {
95 /// Time points
96 pub t: Array1<f64>,
97
98 /// Solution values, indexed as [time, space...]
99 pub u: Vec<Array2<f64>>,
100
101 /// Solver information
102 pub info: Option<String>,
103
104 /// Computation time
105 pub computation_time: f64,
106
107 /// Number of time steps
108 pub num_steps: usize,
109
110 /// Number of linear system solves
111 pub num_linear_solves: usize,
112}
113
114/// Crank-Nicolson solver for 1D parabolic PDEs
115///
116/// This solver uses the Crank-Nicolson method, which is a second-order
117/// implicit method that is unconditionally stable for the heat equation.
118pub struct CrankNicolson1D {
119 /// Spatial domain
120 domain: Domain,
121
122 /// Time range [t_start, t_end]
123 time_range: [f64; 2],
124
125 /// Diffusion coefficient function D(x, t, u) for ∂²u/∂x²
126 diffusion_coeff: CoefficientFunction,
127
128 /// Advection coefficient function v(x, t, u) for ∂u/∂x
129 advection_coeff: Option<CoefficientFunction>,
130
131 /// Reaction term function f(x, t, u)
132 reaction_term: Option<CoefficientFunction>,
133
134 /// Initial condition function u(x, 0)
135 initial_condition: InitialCondition,
136
137 /// Boundary conditions
138 boundary_conditions: Vec<BoundaryCondition<f64>>,
139
140 /// Finite difference scheme for spatial discretization
141 fd_scheme: FiniteDifferenceScheme,
142
143 /// Solver options
144 options: ImplicitOptions,
145}
146
147impl CrankNicolson1D {
148 /// Create a new Crank-Nicolson solver for 1D parabolic PDEs
149 pub fn new(
150 domain: Domain,
151 time_range: [f64; 2],
152 diffusion_coeff: impl Fn(f64, f64, f64) -> f64 + Send + Sync + 'static,
153 initial_condition: impl Fn(f64) -> f64 + Send + Sync + 'static,
154 boundary_conditions: Vec<BoundaryCondition<f64>>,
155 options: Option<ImplicitOptions>,
156 ) -> PDEResult<Self> {
157 // Validate the domain
158 if domain.dimensions() != 1 {
159 return Err(PDEError::DomainError(
160 "Domain must be 1-dimensional for 1D Crank-Nicolson solver".to_string(),
161 ));
162 }
163
164 // Validate time _range
165 if time_range[0] >= time_range[1] {
166 return Err(PDEError::DomainError(
167 "Invalid time _range: start must be less than end".to_string(),
168 ));
169 }
170
171 // Validate boundary _conditions
172 if boundary_conditions.len() != 2 {
173 return Err(PDEError::BoundaryConditions(
174 "1D parabolic PDE requires exactly 2 boundary _conditions".to_string(),
175 ));
176 }
177
178 // Ensure we have both lower and upper boundary _conditions
179 let has_lower = boundary_conditions
180 .iter()
181 .any(|bc| bc.location == BoundaryLocation::Lower);
182 let has_upper = boundary_conditions
183 .iter()
184 .any(|bc| bc.location == BoundaryLocation::Upper);
185
186 if !has_lower || !has_upper {
187 return Err(PDEError::BoundaryConditions(
188 "1D parabolic PDE requires both lower and upper boundary _conditions".to_string(),
189 ));
190 }
191
192 let mut options = options.unwrap_or_default();
193 options.method = ImplicitMethod::CrankNicolson;
194
195 Ok(CrankNicolson1D {
196 domain,
197 time_range,
198 diffusion_coeff: Box::new(diffusion_coeff),
199 advection_coeff: None,
200 reaction_term: None,
201 initial_condition: Box::new(initial_condition),
202 boundary_conditions,
203 fd_scheme: FiniteDifferenceScheme::CentralDifference,
204 options,
205 })
206 }
207
208 /// Add an advection term to the PDE
209 pub fn with_advection(
210 mut self,
211 advection_coeff: impl Fn(f64, f64, f64) -> f64 + Send + Sync + 'static,
212 ) -> Self {
213 self.advection_coeff = Some(Box::new(advection_coeff));
214 self
215 }
216
217 /// Add a reaction term to the PDE
218 pub fn with_reaction(
219 mut self,
220 reaction_term: impl Fn(f64, f64, f64) -> f64 + Send + Sync + 'static,
221 ) -> Self {
222 self.reaction_term = Some(Box::new(reaction_term));
223 self
224 }
225
226 /// Set the finite difference scheme for spatial discretization
227 pub fn with_fd_scheme(mut self, scheme: FiniteDifferenceScheme) -> Self {
228 self.fd_scheme = scheme;
229 self
230 }
231
232 /// Solve the PDE using the Crank-Nicolson method
233 pub fn solve(&self) -> PDEResult<ImplicitResult> {
234 let start_time = Instant::now();
235
236 // Generate spatial grid
237 let x_grid = self.domain.grid(0)?;
238 let nx = x_grid.len();
239 let dx = self.domain.grid_spacing(0)?;
240
241 // Time step
242 let dt = self.options.dt.unwrap_or(0.01);
243
244 // Calculate number of time steps
245 let t_start = self.time_range[0];
246 let t_end = self.time_range[1];
247 let num_steps = ((t_end - t_start) / dt).ceil() as usize;
248
249 // Initialize time array
250 let mut t_values = Vec::with_capacity(num_steps + 1);
251 t_values.push(t_start);
252
253 // Initialize solution arrays
254 let mut u_current = Array1::zeros(nx);
255
256 // Apply initial condition
257 for (i, &x) in x_grid.iter().enumerate() {
258 u_current[i] = (self.initial_condition)(x);
259 }
260
261 // Apply Dirichlet boundary conditions to initial condition
262 apply_dirichlet_conditions_to_initial_1d(&mut u_current, &self.boundary_conditions);
263
264 // Store solutions
265 let save_every = self.options.save_every.unwrap_or(1);
266 let mut solutions = Vec::with_capacity((num_steps + 1) / save_every + 1);
267 solutions.push(u_current.clone().into_shape_with_order((nx, 1)).unwrap());
268
269 // Initialize matrices for Crank-Nicolson method
270 let mut a_matrix = Array2::zeros((nx, nx));
271 let mut b_matrix = Array2::zeros((nx, nx));
272
273 // Track solver statistics
274 let mut num_linear_solves = 0;
275
276 // Time-stepping loop
277 for step in 0..num_steps {
278 let t_current = t_start + step as f64 * dt;
279 let t_next = t_current + dt;
280
281 // Set up coefficient matrices based on PDE terms and boundary conditions
282 self.setup_coefficient_matrices(
283 &mut a_matrix,
284 &mut b_matrix,
285 &x_grid,
286 dx,
287 dt,
288 t_current,
289 );
290
291 // Right-hand side vector
292 let rhs = b_matrix.dot(&u_current);
293
294 // Solve the linear system: A * u_{n+1} = B * u_n
295 let u_next = self.solve_linear_system(&a_matrix, &rhs.view())?;
296 num_linear_solves += 1;
297
298 // Update current solution and time
299 u_current = u_next;
300 t_values.push(t_next);
301
302 // Save solution if needed
303 if (step + 1) % save_every == 0 || step == num_steps - 1 {
304 solutions.push(u_current.clone().into_shape_with_order((nx, 1)).unwrap());
305 }
306
307 // Print progress if verbose
308 if self.options.verbose && (step + 1) % 10 == 0 {
309 println!(
310 "Step {}/{} completed, t = {:.4}",
311 step + 1,
312 num_steps,
313 t_next
314 );
315 }
316 }
317
318 // Convert time values to Array1
319 let t_array = Array1::from_vec(t_values);
320
321 // Compute solution time
322 let computation_time = start_time.elapsed().as_secs_f64();
323
324 // Create result
325 let info = Some(format!(
326 "Time steps: {num_steps}, Linear system solves: {num_linear_solves}"
327 ));
328
329 Ok(ImplicitResult {
330 t: t_array,
331 u: solutions,
332 info,
333 computation_time,
334 num_steps,
335 num_linear_solves,
336 })
337 }
338
339 /// Set up coefficient matrices for the Crank-Nicolson method
340 fn setup_coefficient_matrices(
341 &self,
342 a_matrix: &mut Array2<f64>,
343 b_matrix: &mut Array2<f64>,
344 x_grid: &Array1<f64>,
345 dx: f64,
346 dt: f64,
347 t: f64,
348 ) {
349 let nx = x_grid.len();
350
351 // Clear matrices
352 a_matrix.fill(0.0);
353 b_matrix.fill(0.0);
354
355 // Set up implicit (A) and explicit (B) matrices for interior points
356 for i in 1..nx - 1 {
357 let x = x_grid[i];
358 let u_val = 0.0; // Used for linearization around previous state if needed
359
360 // Diffusion coefficient at the current point
361 let d = (self.diffusion_coeff)(x, t, u_val);
362
363 // Crank-Nicolson coefficients for diffusion term
364 let r = 0.5 * d * dt / (dx * dx);
365
366 // Implicit part (left-hand side)
367 a_matrix[[i, i - 1]] = -r; // Coefficient for u_{i-1}^{n+1}
368 a_matrix[[i, i]] = 1.0 + 2.0 * r; // Coefficient for u_{i}^{n+1}
369 a_matrix[[i, i + 1]] = -r; // Coefficient for u_{i+1}^{n+1}
370
371 // Explicit part (right-hand side)
372 b_matrix[[i, i - 1]] = r; // Coefficient for u_{i-1}^{n}
373 b_matrix[[i, i]] = 1.0 - 2.0 * r; // Coefficient for u_{i}^{n}
374 b_matrix[[i, i + 1]] = r; // Coefficient for u_{i+1}^{n}
375
376 // Add advection term if present
377 if let Some(advection) = &self.advection_coeff {
378 let v = advection(x, t, u_val);
379
380 // Coefficient for advection term (central difference)
381 let c = 0.25 * v * dt / dx;
382
383 // Implicit part
384 a_matrix[[i, i - 1]] -= c; // Additional term for u_{i-1}^{n+1}
385 a_matrix[[i, i + 1]] += c; // Additional term for u_{i+1}^{n+1}
386
387 // Explicit part
388 b_matrix[[i, i - 1]] -= c; // Additional term for u_{i-1}^{n}
389 b_matrix[[i, i + 1]] += c; // Additional term for u_{i+1}^{n}
390 }
391
392 // Add reaction term if present
393 if let Some(reaction) = &self.reaction_term {
394 // For linear reaction terms of the form R(u) = ku
395 // For nonlinear terms, would need to linearize around previous state
396 let k = reaction(x, t, u_val);
397
398 // Coefficient for reaction term
399 let s = 0.5 * k * dt;
400
401 // Implicit part
402 a_matrix[[i, i]] += s; // Additional term for u_{i}^{n+1}
403
404 // Explicit part
405 b_matrix[[i, i]] += s; // Additional term for u_{i}^{n}
406 }
407 }
408
409 // Apply boundary conditions
410 for bc in &self.boundary_conditions {
411 match bc.location {
412 BoundaryLocation::Lower => {
413 // Apply boundary condition at x[0]
414 match bc.bc_type {
415 BoundaryConditionType::Dirichlet => {
416 // u(a, t) = bc.value
417 // Set the row to enforce u_0^{n+1} = bc.value
418 for j in 0..nx {
419 a_matrix[[0, j]] = 0.0;
420 b_matrix[[0, j]] = 0.0;
421 }
422 a_matrix[[0, 0]] = 1.0;
423 b_matrix[[0, 0]] = 0.0;
424
425 // Add boundary value to RHS
426 b_matrix[[0, 0]] = bc.value;
427 }
428 BoundaryConditionType::Neumann => {
429 // du/dx(a, t) = bc.value
430 // Use second-order one-sided difference:
431 // (-3u_0 + 4u_1 - u_2)/(2dx) = bc.value
432
433 for j in 0..nx {
434 a_matrix[[0, j]] = 0.0;
435 b_matrix[[0, j]] = 0.0;
436 }
437
438 // Implicit part
439 a_matrix[[0, 0]] = -3.0;
440 a_matrix[[0, 1]] = 4.0;
441 a_matrix[[0, 2]] = -1.0;
442
443 // RHS is the boundary value scaled by 2*dx
444 b_matrix[[0, 0]] = 2.0 * dx * bc.value;
445 }
446 BoundaryConditionType::Robin => {
447 // a*u + b*du/dx = c
448 if let Some([a_val, b_val, c_val]) = bc.coefficients {
449 // Use second-order one-sided difference for the derivative:
450 // (-3u_0 + 4u_1 - u_2)/(2dx)
451
452 for j in 0..nx {
453 a_matrix[[0, j]] = 0.0;
454 b_matrix[[0, j]] = 0.0;
455 }
456
457 // Implicit part: a*u_0 + b*(-3u_0 + 4u_1 - u_2)/(2dx) = c
458 a_matrix[[0, 0]] = a_val - 3.0 * b_val / (2.0 * dx);
459 a_matrix[[0, 1]] = 4.0 * b_val / (2.0 * dx);
460 a_matrix[[0, 2]] = -b_val / (2.0 * dx);
461
462 // RHS is the boundary value c
463 b_matrix[[0, 0]] = c_val;
464 }
465 }
466 BoundaryConditionType::Periodic => {
467 // For periodic BCs, special handling at both boundaries together is needed
468 // (see below)
469 }
470 }
471 }
472 BoundaryLocation::Upper => {
473 // Apply boundary condition at x[nx-1]
474 let i = nx - 1;
475
476 match bc.bc_type {
477 BoundaryConditionType::Dirichlet => {
478 // u(b, t) = bc.value
479 // Set the row to enforce u_{nx-1}^{n+1} = bc.value
480 for j in 0..nx {
481 a_matrix[[i, j]] = 0.0;
482 b_matrix[[i, j]] = 0.0;
483 }
484 a_matrix[[i, i]] = 1.0;
485 b_matrix[[i, i]] = 0.0;
486
487 // Add boundary value to RHS
488 b_matrix[[i, i]] = bc.value;
489 }
490 BoundaryConditionType::Neumann => {
491 // du/dx(b, t) = bc.value
492 // Use second-order one-sided difference:
493 // (3u_{nx-1} - 4u_{nx-2} + u_{nx-3})/(2dx) = bc.value
494
495 for j in 0..nx {
496 a_matrix[[i, j]] = 0.0;
497 b_matrix[[i, j]] = 0.0;
498 }
499
500 // Implicit part
501 a_matrix[[i, i]] = 3.0;
502 a_matrix[[i, i - 1]] = -4.0;
503 a_matrix[[i, i - 2]] = 1.0;
504
505 // RHS is the boundary value scaled by 2*dx
506 b_matrix[[i, i]] = 2.0 * dx * bc.value;
507 }
508 BoundaryConditionType::Robin => {
509 // a*u + b*du/dx = c
510 if let Some([a_val, b_val, c_val]) = bc.coefficients {
511 // Use second-order one-sided difference for the derivative:
512 // (3u_{nx-1} - 4u_{nx-2} + u_{nx-3})/(2dx)
513
514 for j in 0..nx {
515 a_matrix[[i, j]] = 0.0;
516 b_matrix[[i, j]] = 0.0;
517 }
518
519 // Implicit part: a*u_{nx-1} + b*(3u_{nx-1} - 4u_{nx-2} + u_{nx-3})/(2dx) = c
520 a_matrix[[i, i]] = a_val + 3.0 * b_val / (2.0 * dx);
521 a_matrix[[i, i - 1]] = -4.0 * b_val / (2.0 * dx);
522 a_matrix[[i, i - 2]] = b_val / (2.0 * dx);
523
524 // RHS is the boundary value c
525 b_matrix[[i, i]] = c_val;
526 }
527 }
528 BoundaryConditionType::Periodic => {
529 // Handle periodic boundary conditions
530 // We need to make u_0 = u_{nx-1} and ensure the fluxes match at the boundaries
531
532 // First, clear the boundary rows
533 for j in 0..nx {
534 a_matrix[[0, j]] = 0.0;
535 a_matrix[[i, j]] = 0.0;
536 b_matrix[[0, j]] = 0.0;
537 b_matrix[[i, j]] = 0.0;
538 }
539
540 // For the lower boundary (i=0), set up equation connecting to upper boundary
541 // Use periodic stencil for diffusion
542 let x = x_grid[0];
543 let u_val = 0.0; // Used for linearization if needed
544 let d = (self.diffusion_coeff)(x, t, u_val);
545 let r = 0.5 * d * dt / (dx * dx);
546
547 a_matrix[[0, i]] = -r; // Coefficient for u_{nx-1}^{n+1}
548 a_matrix[[0, 0]] = 1.0 + 2.0 * r; // Coefficient for u_{0}^{n+1}
549 a_matrix[[0, 1]] = -r; // Coefficient for u_{1}^{n+1}
550
551 b_matrix[[0, i]] = r; // Coefficient for u_{nx-1}^{n}
552 b_matrix[[0, 0]] = 1.0 - 2.0 * r; // Coefficient for u_{0}^{n}
553 b_matrix[[0, 1]] = r; // Coefficient for u_{1}^{n}
554
555 // For the upper boundary (i=nx-1), set up equation connecting to lower boundary
556 let x = x_grid[i];
557 let d = (self.diffusion_coeff)(x, t, u_val);
558 let r = 0.5 * d * dt / (dx * dx);
559
560 a_matrix[[i, i - 1]] = -r; // Coefficient for u_{nx-2}^{n+1}
561 a_matrix[[i, i]] = 1.0 + 2.0 * r; // Coefficient for u_{nx-1}^{n+1}
562 a_matrix[[i, 0]] = -r; // Coefficient for u_{0}^{n+1}
563
564 b_matrix[[i, i - 1]] = r; // Coefficient for u_{nx-2}^{n}
565 b_matrix[[i, i]] = 1.0 - 2.0 * r; // Coefficient for u_{nx-1}^{n}
566 b_matrix[[i, 0]] = r; // Coefficient for u_{0}^{n}
567
568 // Add advection term if present
569 if let Some(advection) = &self.advection_coeff {
570 // Lower boundary
571 let v = advection(x_grid[0], t, u_val);
572 let c = 0.25 * v * dt / dx;
573
574 a_matrix[[0, i]] -= c; // Additional term for u_{nx-1}^{n+1}
575 a_matrix[[0, 1]] += c; // Additional term for u_{1}^{n+1}
576
577 b_matrix[[0, i]] -= c; // Additional term for u_{nx-1}^{n}
578 b_matrix[[0, 1]] += c; // Additional term for u_{1}^{n}
579
580 // Upper boundary
581 let v = advection(x_grid[i], t, u_val);
582 let c = 0.25 * v * dt / dx;
583
584 a_matrix[[i, i - 1]] -= c; // Additional term for u_{nx-2}^{n+1}
585 a_matrix[[i, 0]] += c; // Additional term for u_{0}^{n+1}
586
587 b_matrix[[i, i - 1]] -= c; // Additional term for u_{nx-2}^{n}
588 b_matrix[[i, 0]] += c; // Additional term for u_{0}^{n}
589 }
590 }
591 }
592 }
593 }
594 }
595
596 // Special case for periodic boundary conditions
597 // If we have both lower and upper periodic boundary conditions,
598 // we need to make sure they're consistent
599 let has_periodic_lower = self.boundary_conditions.iter().any(|bc| {
600 bc.location == BoundaryLocation::Lower && bc.bc_type == BoundaryConditionType::Periodic
601 });
602
603 let has_periodic_upper = self.boundary_conditions.iter().any(|bc| {
604 bc.location == BoundaryLocation::Upper && bc.bc_type == BoundaryConditionType::Periodic
605 });
606
607 if has_periodic_lower && has_periodic_upper {
608 // Already handled in the boundary condition loop
609 }
610 }
611
612 /// Solve the linear system Ax = b
613 fn solve_linear_system(&self, a: &Array2<f64>, b: &ArrayView1<f64>) -> PDEResult<Array1<f64>> {
614 let n = b.len();
615
616 // Simple tridiagonal solver for Crank-Nicolson matrices
617 // For a general solver, use a specialized linear algebra library
618
619 // Check if the matrix is tridiagonal
620 let is_tridiagonal = a
621 .indexed_iter()
622 .filter(|((i, j), &val)| val != 0.0 && (*i as isize - *j as isize).abs() > 1)
623 .count()
624 == 0;
625
626 if is_tridiagonal {
627 // Extract the tridiagonal elements
628 let mut lower = Array1::zeros(n - 1);
629 let mut diagonal = Array1::zeros(n);
630 let mut upper = Array1::zeros(n - 1);
631
632 for i in 0..n {
633 diagonal[i] = a[[i, i]];
634 if i < n - 1 {
635 upper[i] = a[[i, i + 1]];
636 }
637 if i > 0 {
638 lower[i - 1] = a[[i, i - 1]];
639 }
640 }
641
642 // Solve tridiagonal system using Thomas algorithm
643 let mut solution = Array1::zeros(n);
644 let mut temp_diag = diagonal.clone();
645 let mut temp_rhs = b.to_owned();
646
647 // Forward sweep
648 for i in 1..n {
649 let w = lower[i - 1] / temp_diag[i - 1];
650 temp_diag[i] -= w * upper[i - 1];
651 temp_rhs[i] -= w * temp_rhs[i - 1];
652 }
653
654 // Backward sweep
655 solution[n - 1] = temp_rhs[n - 1] / temp_diag[n - 1];
656 for i in (0..n - 1).rev() {
657 solution[i] = (temp_rhs[i] - upper[i] * solution[i + 1]) / temp_diag[i];
658 }
659
660 Ok(solution)
661 } else {
662 // General case: Gaussian elimination with partial pivoting
663 // For a real implementation, use a specialized linear algebra library
664
665 // Create copies of A and b
666 let mut a_copy = a.clone();
667 let mut b_copy = b.to_owned();
668
669 // Forward elimination
670 for i in 0..n {
671 // Find pivot
672 let mut max_val = a_copy[[i, i]].abs();
673 let mut max_row = i;
674
675 for k in i + 1..n {
676 if a_copy[[k, i]].abs() > max_val {
677 max_val = a_copy[[k, i]].abs();
678 max_row = k;
679 }
680 }
681
682 // Check if matrix is singular
683 if max_val < 1e-10 {
684 return Err(PDEError::Other(
685 "Matrix is singular or nearly singular".to_string(),
686 ));
687 }
688
689 // Swap rows if necessary
690 if max_row != i {
691 for j in i..n {
692 let temp = a_copy[[i, j]];
693 a_copy[[i, j]] = a_copy[[max_row, j]];
694 a_copy[[max_row, j]] = temp;
695 }
696
697 let temp = b_copy[i];
698 b_copy[i] = b_copy[max_row];
699 b_copy[max_row] = temp;
700 }
701
702 // Eliminate below
703 for k in i + 1..n {
704 let factor = a_copy[[k, i]] / a_copy[[i, i]];
705
706 for j in i..n {
707 a_copy[[k, j]] -= factor * a_copy[[i, j]];
708 }
709
710 b_copy[k] -= factor * b_copy[i];
711 }
712 }
713
714 // Back substitution
715 let mut x = Array1::zeros(n);
716 for i in (0..n).rev() {
717 let mut sum = 0.0;
718 for j in i + 1..n {
719 sum += a_copy[[i, j]] * x[j];
720 }
721
722 x[i] = (b_copy[i] - sum) / a_copy[[i, i]];
723 }
724
725 Ok(x)
726 }
727 }
728}
729
730/// Backward Euler solver for 1D parabolic PDEs
731///
732/// This solver uses the Backward Euler method, which is a first-order
733/// fully implicit method that is L-stable and well-suited for stiff problems.
734pub struct BackwardEuler1D {
735 /// Spatial domain
736 domain: Domain,
737
738 /// Time range [t_start, t_end]
739 time_range: [f64; 2],
740
741 /// Diffusion coefficient function D(x, t, u) for ∂²u/∂x²
742 diffusion_coeff: CoefficientFunction,
743
744 /// Advection coefficient function v(x, t, u) for ∂u/∂x
745 advection_coeff: Option<CoefficientFunction>,
746
747 /// Reaction term function f(x, t, u)
748 reaction_term: Option<CoefficientFunction>,
749
750 /// Initial condition function u(x, 0)
751 initial_condition: InitialCondition,
752
753 /// Boundary conditions
754 boundary_conditions: Vec<BoundaryCondition<f64>>,
755
756 /// Finite difference scheme for spatial discretization
757 fd_scheme: FiniteDifferenceScheme,
758
759 /// Solver options
760 options: ImplicitOptions,
761}
762
763impl BackwardEuler1D {
764 /// Create a new Backward Euler solver for 1D parabolic PDEs
765 pub fn new(
766 domain: Domain,
767 time_range: [f64; 2],
768 diffusion_coeff: impl Fn(f64, f64, f64) -> f64 + Send + Sync + 'static,
769 initial_condition: impl Fn(f64) -> f64 + Send + Sync + 'static,
770 boundary_conditions: Vec<BoundaryCondition<f64>>,
771 options: Option<ImplicitOptions>,
772 ) -> PDEResult<Self> {
773 // Validate the domain
774 if domain.dimensions() != 1 {
775 return Err(PDEError::DomainError(
776 "Domain must be 1-dimensional for 1D Backward Euler solver".to_string(),
777 ));
778 }
779
780 // Validate time _range
781 if time_range[0] >= time_range[1] {
782 return Err(PDEError::DomainError(
783 "Invalid time _range: start must be less than end".to_string(),
784 ));
785 }
786
787 // Validate boundary _conditions
788 if boundary_conditions.len() != 2 {
789 return Err(PDEError::BoundaryConditions(
790 "1D parabolic PDE requires exactly 2 boundary _conditions".to_string(),
791 ));
792 }
793
794 // Ensure we have both lower and upper boundary _conditions
795 let has_lower = boundary_conditions
796 .iter()
797 .any(|bc| bc.location == BoundaryLocation::Lower);
798 let has_upper = boundary_conditions
799 .iter()
800 .any(|bc| bc.location == BoundaryLocation::Upper);
801
802 if !has_lower || !has_upper {
803 return Err(PDEError::BoundaryConditions(
804 "1D parabolic PDE requires both lower and upper boundary _conditions".to_string(),
805 ));
806 }
807
808 let mut options = options.unwrap_or_default();
809 options.method = ImplicitMethod::BackwardEuler;
810
811 Ok(BackwardEuler1D {
812 domain,
813 time_range,
814 diffusion_coeff: Box::new(diffusion_coeff),
815 advection_coeff: None,
816 reaction_term: None,
817 initial_condition: Box::new(initial_condition),
818 boundary_conditions,
819 fd_scheme: FiniteDifferenceScheme::CentralDifference,
820 options,
821 })
822 }
823
824 /// Add an advection term to the PDE
825 pub fn with_advection(
826 mut self,
827 advection_coeff: impl Fn(f64, f64, f64) -> f64 + Send + Sync + 'static,
828 ) -> Self {
829 self.advection_coeff = Some(Box::new(advection_coeff));
830 self
831 }
832
833 /// Add a reaction term to the PDE
834 pub fn with_reaction(
835 mut self,
836 reaction_term: impl Fn(f64, f64, f64) -> f64 + Send + Sync + 'static,
837 ) -> Self {
838 self.reaction_term = Some(Box::new(reaction_term));
839 self
840 }
841
842 /// Set the finite difference scheme for spatial discretization
843 pub fn with_fd_scheme(mut self, scheme: FiniteDifferenceScheme) -> Self {
844 self.fd_scheme = scheme;
845 self
846 }
847
848 /// Solve the PDE using the Backward Euler method
849 pub fn solve(&self) -> PDEResult<ImplicitResult> {
850 let start_time = Instant::now();
851
852 // Generate spatial grid
853 let x_grid = self.domain.grid(0)?;
854 let nx = x_grid.len();
855 let dx = self.domain.grid_spacing(0)?;
856
857 // Time step
858 let dt = self.options.dt.unwrap_or(0.01);
859
860 // Calculate number of time steps
861 let t_start = self.time_range[0];
862 let t_end = self.time_range[1];
863 let num_steps = ((t_end - t_start) / dt).ceil() as usize;
864
865 // Initialize time array
866 let mut t_values = Vec::with_capacity(num_steps + 1);
867 t_values.push(t_start);
868
869 // Initialize solution arrays
870 let mut u_current = Array1::zeros(nx);
871
872 // Apply initial condition
873 for (i, &x) in x_grid.iter().enumerate() {
874 u_current[i] = (self.initial_condition)(x);
875 }
876
877 // Apply Dirichlet boundary conditions to initial condition
878 apply_dirichlet_conditions_to_initial_1d(&mut u_current, &self.boundary_conditions);
879
880 // Store solutions
881 let save_every = self.options.save_every.unwrap_or(1);
882 let mut solutions = Vec::with_capacity((num_steps + 1) / save_every + 1);
883 solutions.push(u_current.clone().into_shape_with_order((nx, 1)).unwrap());
884
885 // Initialize coefficient matrix for Backward Euler method
886 let mut a_matrix = Array2::zeros((nx, nx));
887
888 // Track solver statistics
889 let mut num_linear_solves = 0;
890
891 // Time-stepping loop
892 for step in 0..num_steps {
893 let t_current = t_start + step as f64 * dt;
894 let t_next = t_current + dt;
895
896 // Set up coefficient matrix based on PDE terms and boundary conditions
897 self.setup_coefficient_matrix(
898 &mut a_matrix,
899 &x_grid,
900 dx,
901 dt,
902 t_next, // Use t_{n+1} for fully implicit scheme
903 );
904
905 // Right-hand side vector is just the current solution
906 let rhs = u_current.clone();
907
908 // Solve the linear system: A * u_{n+1} = u_n
909 let u_next = self.solve_linear_system(&a_matrix, &rhs.view())?;
910 num_linear_solves += 1;
911
912 // Update current solution and time
913 u_current = u_next;
914 t_values.push(t_next);
915
916 // Save solution if needed
917 if (step + 1) % save_every == 0 || step == num_steps - 1 {
918 solutions.push(u_current.clone().into_shape_with_order((nx, 1)).unwrap());
919 }
920
921 // Print progress if verbose
922 if self.options.verbose && (step + 1) % 10 == 0 {
923 println!(
924 "Step {}/{} completed, t = {:.4}",
925 step + 1,
926 num_steps,
927 t_next
928 );
929 }
930 }
931
932 // Convert time values to Array1
933 let t_array = Array1::from_vec(t_values);
934
935 // Compute solution time
936 let computation_time = start_time.elapsed().as_secs_f64();
937
938 // Create result
939 let info = Some(format!(
940 "Time steps: {num_steps}, Linear system solves: {num_linear_solves}"
941 ));
942
943 Ok(ImplicitResult {
944 t: t_array,
945 u: solutions,
946 info,
947 computation_time,
948 num_steps,
949 num_linear_solves,
950 })
951 }
952
953 /// Set up coefficient matrix for the Backward Euler method
954 fn setup_coefficient_matrix(
955 &self,
956 a_matrix: &mut Array2<f64>,
957 x_grid: &Array1<f64>,
958 dx: f64,
959 dt: f64,
960 t: f64,
961 ) {
962 let nx = x_grid.len();
963
964 // Clear _matrix
965 a_matrix.fill(0.0);
966
967 // Set up implicit _matrix for interior points
968 for i in 1..nx - 1 {
969 let x = x_grid[i];
970 let u_val = 0.0; // Used for linearization around previous state if needed
971
972 // Diffusion coefficient at the current point
973 let d = (self.diffusion_coeff)(x, t, u_val);
974
975 // Backward Euler coefficients for diffusion term
976 let r = d * dt / (dx * dx);
977
978 // Coefficient _matrix for implicit scheme
979 a_matrix[[i, i - 1]] = -r; // Coefficient for u_{i-1}^{n+1}
980 a_matrix[[i, i]] = 1.0 + 2.0 * r; // Coefficient for u_{i}^{n+1}
981 a_matrix[[i, i + 1]] = -r; // Coefficient for u_{i+1}^{n+1}
982
983 // Add advection term if present
984 if let Some(advection) = &self.advection_coeff {
985 let v = advection(x, t, u_val);
986
987 // Coefficient for advection term (upwind for stability)
988 if v > 0.0 {
989 // Use backward difference for positive velocity
990 let c = v * dt / dx;
991 a_matrix[[i, i]] += c;
992 a_matrix[[i, i - 1]] -= c;
993 } else {
994 // Use forward difference for negative velocity
995 let c = -v * dt / dx;
996 a_matrix[[i, i]] += c;
997 a_matrix[[i, i + 1]] -= c;
998 }
999 }
1000
1001 // Add reaction term if present
1002 if let Some(reaction) = &self.reaction_term {
1003 // For linear reaction terms of the form R(u) = ku
1004 // For nonlinear terms, would need to linearize around previous state
1005 let k = reaction(x, t, u_val);
1006
1007 // Coefficient for reaction term
1008 let s = k * dt;
1009 a_matrix[[i, i]] += s; // Additional term for u_{i}^{n+1}
1010 }
1011 }
1012
1013 // Apply boundary conditions
1014 for bc in &self.boundary_conditions {
1015 match bc.location {
1016 BoundaryLocation::Lower => {
1017 // Apply boundary condition at x[0]
1018 match bc.bc_type {
1019 BoundaryConditionType::Dirichlet => {
1020 // u(a, t) = bc.value
1021 // Replace the equation with u_0^{n+1} = bc.value
1022 for j in 0..nx {
1023 a_matrix[[0, j]] = 0.0;
1024 }
1025 a_matrix[[0, 0]] = 1.0;
1026
1027 // Right-hand side will have the boundary value
1028 // For Backward Euler, we directly modify the solution vector
1029 // to enforce the boundary condition
1030 }
1031 BoundaryConditionType::Neumann => {
1032 // du/dx(a, t) = bc.value
1033 // Use second-order one-sided difference:
1034 // (-3u_0 + 4u_1 - u_2)/(2dx) = bc.value
1035
1036 for j in 0..nx {
1037 a_matrix[[0, j]] = 0.0;
1038 }
1039
1040 a_matrix[[0, 0]] = -3.0;
1041 a_matrix[[0, 1]] = 4.0;
1042 a_matrix[[0, 2]] = -1.0;
1043
1044 // Right-hand side will have the boundary value scaled by 2*dx
1045 // For Backward Euler, we directly modify the solution vector
1046 }
1047 BoundaryConditionType::Robin => {
1048 // a*u + b*du/dx = c
1049 if let Some([a_val, b_val, c_val]) = bc.coefficients {
1050 // Use second-order one-sided difference for the derivative:
1051 // (-3u_0 + 4u_1 - u_2)/(2dx)
1052
1053 for j in 0..nx {
1054 a_matrix[[0, j]] = 0.0;
1055 }
1056
1057 a_matrix[[0, 0]] = a_val - 3.0 * b_val / (2.0 * dx);
1058 a_matrix[[0, 1]] = 4.0 * b_val / (2.0 * dx);
1059 a_matrix[[0, 2]] = -b_val / (2.0 * dx);
1060
1061 // Right-hand side will have the boundary value c
1062 // For Backward Euler, we directly modify the solution vector
1063 }
1064 }
1065 BoundaryConditionType::Periodic => {
1066 // For periodic BCs, special handling at both boundaries together is needed
1067 // (see below)
1068 }
1069 }
1070 }
1071 BoundaryLocation::Upper => {
1072 // Apply boundary condition at x[nx-1]
1073 let i = nx - 1;
1074
1075 match bc.bc_type {
1076 BoundaryConditionType::Dirichlet => {
1077 // u(b, t) = bc.value
1078 // Replace the equation with u_{nx-1}^{n+1} = bc.value
1079 for j in 0..nx {
1080 a_matrix[[i, j]] = 0.0;
1081 }
1082 a_matrix[[i, i]] = 1.0;
1083
1084 // Right-hand side will have the boundary value
1085 // For Backward Euler, we directly modify the solution vector
1086 }
1087 BoundaryConditionType::Neumann => {
1088 // du/dx(b, t) = bc.value
1089 // Use second-order one-sided difference:
1090 // (3u_{nx-1} - 4u_{nx-2} + u_{nx-3})/(2dx) = bc.value
1091
1092 for j in 0..nx {
1093 a_matrix[[i, j]] = 0.0;
1094 }
1095
1096 a_matrix[[i, i]] = 3.0;
1097 a_matrix[[i, i - 1]] = -4.0;
1098 a_matrix[[i, i - 2]] = 1.0;
1099
1100 // Right-hand side will have the boundary value scaled by 2*dx
1101 // For Backward Euler, we directly modify the solution vector
1102 }
1103 BoundaryConditionType::Robin => {
1104 // a*u + b*du/dx = c
1105 if let Some([a_val, b_val, c_val]) = bc.coefficients {
1106 // Use second-order one-sided difference for the derivative:
1107 // (3u_{nx-1} - 4u_{nx-2} + u_{nx-3})/(2dx)
1108
1109 for j in 0..nx {
1110 a_matrix[[i, j]] = 0.0;
1111 }
1112
1113 a_matrix[[i, i]] = a_val + 3.0 * b_val / (2.0 * dx);
1114 a_matrix[[i, i - 1]] = -4.0 * b_val / (2.0 * dx);
1115 a_matrix[[i, i - 2]] = b_val / (2.0 * dx);
1116
1117 // Right-hand side will have the boundary value c
1118 // For Backward Euler, we directly modify the solution vector
1119 }
1120 }
1121 BoundaryConditionType::Periodic => {
1122 // Handle periodic boundary conditions
1123 // We need to make u_0 = u_{nx-1} and ensure the fluxes match at the boundaries
1124
1125 // First, clear the boundary rows
1126 for j in 0..nx {
1127 a_matrix[[0, j]] = 0.0;
1128 a_matrix[[i, j]] = 0.0;
1129 }
1130
1131 // For the lower boundary (i=0), set up equation connecting to upper boundary
1132 // Use periodic stencil for diffusion
1133 let x = x_grid[0];
1134 let u_val = 0.0; // Used for linearization if needed
1135 let d = (self.diffusion_coeff)(x, t, u_val);
1136 let r = d * dt / (dx * dx);
1137
1138 a_matrix[[0, i]] = -r; // Coefficient for u_{nx-1}^{n+1}
1139 a_matrix[[0, 0]] = 1.0 + 2.0 * r; // Coefficient for u_{0}^{n+1}
1140 a_matrix[[0, 1]] = -r; // Coefficient for u_{1}^{n+1}
1141
1142 // For the upper boundary (i=nx-1), set up equation connecting to lower boundary
1143 let x = x_grid[i];
1144 let d = (self.diffusion_coeff)(x, t, u_val);
1145 let r = d * dt / (dx * dx);
1146
1147 a_matrix[[i, i - 1]] = -r; // Coefficient for u_{nx-2}^{n+1}
1148 a_matrix[[i, i]] = 1.0 + 2.0 * r; // Coefficient for u_{nx-1}^{n+1}
1149 a_matrix[[i, 0]] = -r; // Coefficient for u_{0}^{n+1}
1150
1151 // Add advection term if present
1152 if let Some(advection) = &self.advection_coeff {
1153 // Lower boundary
1154 let v = advection(x_grid[0], t, u_val);
1155
1156 if v > 0.0 {
1157 // Use backward difference for positive velocity
1158 let c = v * dt / dx;
1159 a_matrix[[0, 0]] += c;
1160 a_matrix[[0, i]] -= c;
1161 } else {
1162 // Use forward difference for negative velocity
1163 let c = -v * dt / dx;
1164 a_matrix[[0, 0]] += c;
1165 a_matrix[[0, 1]] -= c;
1166 }
1167
1168 // Upper boundary
1169 let v = advection(x_grid[i], t, u_val);
1170
1171 if v > 0.0 {
1172 // Use backward difference for positive velocity
1173 let c = v * dt / dx;
1174 a_matrix[[i, i]] += c;
1175 a_matrix[[i, i - 1]] -= c;
1176 } else {
1177 // Use forward difference for negative velocity
1178 let c = -v * dt / dx;
1179 a_matrix[[i, i]] += c;
1180 a_matrix[[i, 0]] -= c;
1181 }
1182 }
1183 }
1184 }
1185 }
1186 }
1187 }
1188
1189 // Special handling for Dirichlet boundary conditions
1190 // For Backward Euler, we need to modify the RHS vector with the boundary values
1191 // This is done in the solve method when applying the boundary conditions to the solution
1192 }
1193
1194 /// Solve the linear system Ax = b
1195 fn solve_linear_system(&self, a: &Array2<f64>, b: &ArrayView1<f64>) -> PDEResult<Array1<f64>> {
1196 let n = b.len();
1197
1198 // Simple tridiagonal solver for Backward Euler matrices
1199 // For a general solver, use a specialized linear algebra library
1200
1201 // Check if the matrix is tridiagonal
1202 let is_tridiagonal = a
1203 .indexed_iter()
1204 .filter(|((i, j), &val)| val != 0.0 && (*i as isize - *j as isize).abs() > 1)
1205 .count()
1206 == 0;
1207
1208 if is_tridiagonal {
1209 // Extract the tridiagonal elements
1210 let mut lower = Array1::zeros(n - 1);
1211 let mut diagonal = Array1::zeros(n);
1212 let mut upper = Array1::zeros(n - 1);
1213
1214 for i in 0..n {
1215 diagonal[i] = a[[i, i]];
1216 if i < n - 1 {
1217 upper[i] = a[[i, i + 1]];
1218 }
1219 if i > 0 {
1220 lower[i - 1] = a[[i, i - 1]];
1221 }
1222 }
1223
1224 // Solve tridiagonal system using Thomas algorithm
1225 let mut solution = Array1::zeros(n);
1226 let mut temp_diag = diagonal.clone();
1227 let mut temp_rhs = b.to_owned();
1228
1229 // Forward sweep
1230 for i in 1..n {
1231 let w = lower[i - 1] / temp_diag[i - 1];
1232 temp_diag[i] -= w * upper[i - 1];
1233 temp_rhs[i] -= w * temp_rhs[i - 1];
1234 }
1235
1236 // Backward sweep
1237 solution[n - 1] = temp_rhs[n - 1] / temp_diag[n - 1];
1238 for i in (0..n - 1).rev() {
1239 solution[i] = (temp_rhs[i] - upper[i] * solution[i + 1]) / temp_diag[i];
1240 }
1241
1242 // Apply Dirichlet boundary conditions directly
1243 for bc in &self.boundary_conditions {
1244 if bc.bc_type == BoundaryConditionType::Dirichlet {
1245 match bc.location {
1246 BoundaryLocation::Lower => solution[0] = bc.value,
1247 BoundaryLocation::Upper => solution[n - 1] = bc.value,
1248 }
1249 }
1250 }
1251
1252 Ok(solution)
1253 } else {
1254 // General case: Gaussian elimination with partial pivoting
1255 // For a real implementation, use a specialized linear algebra library
1256
1257 // Create copies of A and b
1258 let mut a_copy = a.clone();
1259 let mut b_copy = b.to_owned();
1260
1261 // Forward elimination
1262 for i in 0..n {
1263 // Find pivot
1264 let mut max_val = a_copy[[i, i]].abs();
1265 let mut max_row = i;
1266
1267 for k in i + 1..n {
1268 if a_copy[[k, i]].abs() > max_val {
1269 max_val = a_copy[[k, i]].abs();
1270 max_row = k;
1271 }
1272 }
1273
1274 // Check if matrix is singular
1275 if max_val < 1e-10 {
1276 return Err(PDEError::Other(
1277 "Matrix is singular or nearly singular".to_string(),
1278 ));
1279 }
1280
1281 // Swap rows if necessary
1282 if max_row != i {
1283 for j in i..n {
1284 let temp = a_copy[[i, j]];
1285 a_copy[[i, j]] = a_copy[[max_row, j]];
1286 a_copy[[max_row, j]] = temp;
1287 }
1288
1289 let temp = b_copy[i];
1290 b_copy[i] = b_copy[max_row];
1291 b_copy[max_row] = temp;
1292 }
1293
1294 // Eliminate below
1295 for k in i + 1..n {
1296 let factor = a_copy[[k, i]] / a_copy[[i, i]];
1297
1298 for j in i..n {
1299 a_copy[[k, j]] -= factor * a_copy[[i, j]];
1300 }
1301
1302 b_copy[k] -= factor * b_copy[i];
1303 }
1304 }
1305
1306 // Back substitution
1307 let mut x = Array1::zeros(n);
1308 for i in (0..n).rev() {
1309 let mut sum = 0.0;
1310 for j in i + 1..n {
1311 sum += a_copy[[i, j]] * x[j];
1312 }
1313
1314 x[i] = (b_copy[i] - sum) / a_copy[[i, i]];
1315 }
1316
1317 // Apply Dirichlet boundary conditions directly
1318 for bc in &self.boundary_conditions {
1319 if bc.bc_type == BoundaryConditionType::Dirichlet {
1320 match bc.location {
1321 BoundaryLocation::Lower => x[0] = bc.value,
1322 BoundaryLocation::Upper => x[n - 1] = bc.value,
1323 }
1324 }
1325 }
1326
1327 Ok(x)
1328 }
1329 }
1330}
1331
1332/// Helper function to apply Dirichlet boundary conditions to initial condition
1333#[allow(dead_code)]
1334fn apply_dirichlet_conditions_to_initial_1d(
1335 u0: &mut Array1<f64>,
1336 boundary_conditions: &[BoundaryCondition<f64>],
1337) {
1338 let nx = u0.len();
1339
1340 for bc in boundary_conditions {
1341 if bc.bc_type == BoundaryConditionType::Dirichlet {
1342 match bc.location {
1343 BoundaryLocation::Lower => {
1344 // Lower boundary (x = x[0])
1345 u0[0] = bc.value;
1346 }
1347 BoundaryLocation::Upper => {
1348 // Upper boundary (x = x[nx-1])
1349 u0[nx - 1] = bc.value;
1350 }
1351 }
1352 }
1353 }
1354}
1355
1356/// Convert an ImplicitResult to a PDESolution
1357impl From<ImplicitResult> for PDESolution<f64> {
1358 fn from(result: ImplicitResult) -> Self {
1359 let mut grids = Vec::new();
1360
1361 // Add time grid
1362 grids.push(result.t.clone());
1363
1364 // Extract spatial grid from solution shape
1365 let nx = result.u[0].shape()[0];
1366
1367 // Create spatial grid (we don't have the actual grid values, so use linspace)
1368 let x_grid = Array1::linspace(0.0, 1.0, nx);
1369 grids.push(x_grid);
1370
1371 // Create solver info
1372 let info = PDESolverInfo {
1373 num_iterations: result.num_linear_solves,
1374 computation_time: result.computation_time,
1375 residual_norm: None,
1376 convergence_history: None,
1377 method: "Implicit Method".to_string(),
1378 };
1379
1380 PDESolution {
1381 grids,
1382 values: result.u,
1383 error_estimate: None,
1384 info,
1385 }
1386 }
1387}