differentiation-schemes

Original source / Raw SKILL.md

---
name: differentiation-schemes
description: Select and apply numerical differentiation schemes for PDE/ODE discretization. Use when choosing finite difference/volume/spectral schemes, building stencils, handling boundaries, estimating truncation error, or analyzing dispersion and dissipation.
allowed-tools: Read, Bash, Write, Grep, Glob
---

# Differentiation Schemes

## Goal

Provide a reliable workflow to select a differentiation scheme, generate stencils, and assess accuracy for simulation discretization.

## Requirements

- Python 3.8+
- NumPy (for stencil computations)
- No heavy dependencies

## Inputs to Gather

| Input | Description | Example |
|-------|-------------|---------|
| Derivative order | First, second, etc. | `1` or `2` |
| Target accuracy | Order of truncation error | `2` or `4` |
| Grid type | Uniform, nonuniform | `uniform` |
| Boundary type | Periodic, Dirichlet, Neumann | `periodic` |
| Smoothness | Smooth or discontinuous | `smooth` |

## Decision Guidance

### Scheme Selection Flowchart

```
Is the field smooth?
├── YES → Is domain periodic?
│   ├── YES → Use central differences or spectral
│   └── NO → Use central interior + one-sided at boundaries
└── NO → Are there shocks/discontinuities?
    ├── YES → Use upwind, TVD, or WENO
    └── NO → Use central with limiters
```

### Quick Reference

| Situation | Recommended Scheme |
|-----------|-------------------|
| Smooth, periodic | Central, spectral |
| Smooth, bounded | Central + one-sided BCs |
| Advection-dominated | Upwind |
| Shocks/fronts | TVD, WENO |
| High accuracy needed | Compact (Padé), spectral |

## Script Outputs (JSON Fields)

| Script | Key Outputs |
|--------|-------------|
| `scripts/stencil_generator.py` | `offsets`, `coefficients`, `order`, `accuracy` |
| `scripts/scheme_selector.py` | `recommended`, `alternatives`, `notes` |
| `scripts/truncation_error.py` | `error_scale`, `order`, `notes` |

## Workflow

1. **Identify requirements** - derivative order, accuracy, smoothness
2. **Select scheme** - Run `scripts/scheme_selector.py`
3. **Generate stencils** - Run `scripts/stencil_generator.py`
4. **Estimate error** - Run `scripts/truncation_error.py`
5. **Validate** - Test with manufactured solutions or grid refinement

## Conversational Workflow Example

**User**: I need to discretize a second derivative for a diffusion equation on a uniform grid. I want 4th-order accuracy.

**Agent workflow**:
1. Select appropriate scheme:
   ```bash
   python3 scripts/scheme_selector.py --smooth --periodic --order 2 --accuracy 4 --json
   ```
2. Generate the stencil:
   ```bash
   python3 scripts/stencil_generator.py --order 2 --accuracy 4 --scheme central --json
   ```
3. Result: 5-point stencil with coefficients `[-1/12, 4/3, -5/2, 4/3, -1/12]` / dx².

## Pre-Discretization Checklist

- [ ] Confirm derivative order and target accuracy
- [ ] Choose scheme appropriate to smoothness and boundaries
- [ ] Generate and inspect stencils at boundaries
- [ ] Estimate truncation error vs physics scales
- [ ] Verify with grid refinement study

## CLI Examples

```bash
# Select scheme for smooth periodic problem
python3 scripts/scheme_selector.py --smooth --periodic --order 1 --accuracy 4 --json

# Generate central difference stencil for first derivative
python3 scripts/stencil_generator.py --order 1 --accuracy 2 --scheme central --json

# Generate 4th-order second derivative stencil
python3 scripts/stencil_generator.py --order 2 --accuracy 4 --scheme central --json

# Estimate truncation error
python3 scripts/truncation_error.py --dx 0.01 --order 2 --accuracy 2 --scale 1.0 --json
```

## Error Handling

| Error | Cause | Resolution |
|-------|-------|------------|
| `order must be positive` | Invalid derivative order | Use 1, 2, 3, ... |
| `accuracy must be even for central` | Odd accuracy requested | Use 2, 4, 6, ... |
| `Unknown scheme` | Invalid scheme type | Use central, upwind, compact |

## Interpretation Guidance

### Stencil Properties

| Property | Meaning |
|----------|---------|
| Symmetric offsets | Central scheme (no directional bias) |
| Asymmetric offsets | One-sided or upwind scheme |
| More points | Higher accuracy but wider stencil |

### Truncation Error Scaling

| Accuracy Order | Error Scales As | Refinement Factor |
|----------------|-----------------|-------------------|
| 2nd order | O(dx²) | 2× refinement → 4× error reduction |
| 4th order | O(dx⁴) | 2× refinement → 16× error reduction |
| 6th order | O(dx⁶) | 2× refinement → 64× error reduction |

### Common Stencils

| Derivative | Accuracy | Points | Coefficients (× 1/dx or 1/dx²) |
|------------|----------|--------|-------------------------------|
| 1st | 2 | 3 | [-1/2, 0, 1/2] |
| 1st | 4 | 5 | [1/12, -2/3, 0, 2/3, -1/12] |
| 2nd | 2 | 3 | [1, -2, 1] |
| 2nd | 4 | 5 | [-1/12, 4/3, -5/2, 4/3, -1/12] |

## Limitations

- **Boundary handling**: Stencil generator provides interior stencils; boundaries need special treatment
- **Nonuniform grids**: Standard stencils assume uniform spacing
- **Spectral**: Not covered by stencil generator

## References

- `references/stencil_catalog.md` - Common stencils
- `references/boundary_handling.md` - One-sided schemes
- `references/scheme_selection.md` - FD/FV/spectral comparison
- `references/error_guidance.md` - Truncation error scaling

## Version History

- **v1.1.0** (2024-12-24): Enhanced documentation, decision guidance, examples
- **v1.0.0**: Initial release with 3 differentiation scripts


---

## Referenced Files

> The following files are referenced in this skill and included for context.

### scripts/scheme_selector.py

```python
#!/usr/bin/env python3
import argparse
import json
import sys
from typing import Dict, List


def select_scheme(
    smooth: bool,
    periodic: bool,
    discontinuous: bool,
    order: int,
    accuracy: int,
    boundary: bool,
) -> Dict[str, object]:
    if order <= 0:
        raise ValueError("order must be positive")
    if accuracy <= 0:
        raise ValueError("accuracy must be positive")
    if discontinuous and smooth:
        raise ValueError("smooth and discontinuous cannot both be true")

    recommended: List[str] = []
    alternatives: List[str] = []
    notes: List[str] = []

    if discontinuous:
        recommended.append("Finite Volume (FV) with limiter/WENO")
        alternatives.append("Upwind FD (low order)")
        notes.append("Avoid high-order central FD near discontinuities.")
    else:
        if periodic and smooth:
            recommended.append("Spectral (Fourier)")
            alternatives.append("High-order central FD")
        elif smooth:
            recommended.append("Central FD")
            if accuracy >= 4:
                alternatives.append("Compact FD")
        else:
            recommended.append("Upwind FD")
            alternatives.append("FV")

    if boundary:
        notes.append("Use one-sided or ghost-cell stencils at boundaries.")

    return {
        "recommended": recommended,
        "alternatives": alternatives,
        "notes": notes,
    }


def parse_args() -> argparse.Namespace:
    parser = argparse.ArgumentParser(
        description="Select a differentiation scheme based on problem characteristics.",
        formatter_class=argparse.ArgumentDefaultsHelpFormatter,
    )
    parser.add_argument("--smooth", action="store_true", help="Field is smooth")
    parser.add_argument("--discontinuous", action="store_true", help="Field has discontinuities")
    parser.add_argument("--periodic", action="store_true", help="Domain is periodic")
    parser.add_argument("--boundary", action="store_true", help="Bounded domain (non-periodic)")
    parser.add_argument("--order", type=int, required=True, help="Derivative order")
    parser.add_argument("--accuracy", type=int, default=2, help="Desired accuracy")
    parser.add_argument("--json", action="store_true", help="Emit JSON output")
    return parser.parse_args()


def main() -> None:
    args = parse_args()
    try:
        result = select_scheme(
            smooth=args.smooth,
            periodic=args.periodic,
            discontinuous=args.discontinuous,
            order=args.order,
            accuracy=args.accuracy,
            boundary=args.boundary,
        )
    except ValueError as exc:
        print(str(exc), file=sys.stderr)
        sys.exit(2)

    payload = {
        "inputs": {
            "smooth": args.smooth,
            "discontinuous": args.discontinuous,
            "periodic": args.periodic,
            "boundary": args.boundary,
            "order": args.order,
            "accuracy": args.accuracy,
        },
        "results": result,
    }

    if args.json:
        print(json.dumps(payload, indent=2, sort_keys=True))
        return

    print("Scheme selection")
    print(f"  recommended: {', '.join(result['recommended'])}")
    if result["alternatives"]:
        print(f"  alternatives: {', '.join(result['alternatives'])}")
    for note in result["notes"]:
        print(f"  note: {note}")


if __name__ == "__main__":
    main()

```

### scripts/stencil_generator.py

```python
#!/usr/bin/env python3
import argparse
import json
import math
import sys
from typing import Dict, List, Optional


def fornberg_coefficients(x: List[float], x0: float, m: int) -> List[float]:
    n = len(x)
    if n == 0:
        raise ValueError("x list must be non-empty")
    if m < 0:
        raise ValueError("m must be non-negative")
    c = [[0.0 for _ in range(m + 1)] for _ in range(n)]
    c[0][0] = 1.0
    c1 = 1.0
    c4 = x[0] - x0
    for i in range(1, n):
        mn = min(i, m)
        c2 = 1.0
        c5 = c4
        c4 = x[i] - x0
        for j in range(i):
            c3 = x[i] - x[j]
            if c3 == 0.0:
                raise ValueError("x points must be distinct")
            c2 *= c3
            if j == i - 1:
                for k in range(mn, 0, -1):
                    c[i][k] = (c1 * (k * c[i - 1][k - 1] - c5 * c[i - 1][k])) / c2
                c[i][0] = (-c1 * c5 * c[i - 1][0]) / c2
            for k in range(mn, 0, -1):
                c[j][k] = (c4 * c[j][k] - k * c[j][k - 1]) / c3
            c[j][0] = (c4 * c[j][0]) / c3
        c1 = c2
    return [c[i][m] for i in range(n)]


def stencil_offsets(order: int, accuracy: int, scheme: str) -> List[int]:
    if order <= 0:
        raise ValueError("order must be positive")
    if accuracy <= 0:
        raise ValueError("accuracy must be positive")
    if scheme not in {"central", "forward", "backward"}:
        raise ValueError("scheme must be central, forward, or backward")

    if scheme == "central":
        points = order + accuracy if order % 2 == 1 else order + accuracy - 1
        if points % 2 == 0:
            points += 1
        radius = points // 2
        return list(range(-radius, radius + 1))

    points = order + accuracy
    if scheme == "forward":
        return list(range(0, points))
    return list(range(0, -points, -1))


def generate_stencil(
    order: int,
    accuracy: int,
    scheme: str,
    dx: float,
    offsets: Optional[List[int]],
) -> Dict[str, object]:
    if dx <= 0:
        raise ValueError("dx must be positive")
    if offsets is None:
        offsets = stencil_offsets(order, accuracy, scheme)
    if not offsets:
        raise ValueError("offsets must be non-empty")

    x = [float(o) * dx for o in offsets]
    coeffs = fornberg_coefficients(x, 0.0, order)
    return {
        "offsets": offsets,
        "coefficients": coeffs,
        "order": order,
        "accuracy": accuracy,
        "scheme": scheme,
    }


def parse_offsets(raw: str) -> List[int]:
    parts = [p.strip() for p in raw.split(",") if p.strip()]
    if not parts:
        raise ValueError("offset list must be a comma-separated list")
    return [int(p) for p in parts]


def parse_args() -> argparse.Namespace:
    parser = argparse.ArgumentParser(
        description="Generate finite difference stencil coefficients.",
        formatter_class=argparse.ArgumentDefaultsHelpFormatter,
    )
    parser.add_argument("--order", type=int, required=True, help="Derivative order")
    parser.add_argument("--accuracy", type=int, default=2, help="Accuracy order")
    parser.add_argument(
        "--scheme",
        choices=["central", "forward", "backward"],
        default="central",
        help="Stencil scheme",
    )
    parser.add_argument("--dx", type=float, default=1.0, help="Grid spacing")
    parser.add_argument("--offsets", default=None, help="Custom offsets (comma-separated)")
    parser.add_argument("--json", action="store_true", help="Emit JSON output")
    return parser.parse_args()


def main() -> None:
    args = parse_args()
    try:
        offsets = parse_offsets(args.offsets) if args.offsets is not None else None
        result = generate_stencil(
            order=args.order,
            accuracy=args.accuracy,
            scheme=args.scheme,
            dx=args.dx,
            offsets=offsets,
        )
    except ValueError as exc:
        print(str(exc), file=sys.stderr)
        sys.exit(2)

    payload = {
        "inputs": {
            "order": args.order,
            "accuracy": args.accuracy,
            "scheme": args.scheme,
            "dx": args.dx,
            "offsets": offsets,
        },
        "results": result,
    }

    if args.json:
        print(json.dumps(payload, indent=2, sort_keys=True))
        return

    print("Stencil")
    print(f"  offsets: {result['offsets']}")
    print(f"  coefficients: {result['coefficients']}")


if __name__ == "__main__":
    main()

```

### scripts/truncation_error.py

```python
#!/usr/bin/env python3
import argparse
import json
import math
import sys
from typing import Dict


def estimate_truncation_error(dx: float, accuracy: int, scale: float) -> Dict[str, object]:
    if dx <= 0:
        raise ValueError("dx must be positive")
    if accuracy <= 0:
        raise ValueError("accuracy must be positive")
    if scale < 0:
        raise ValueError("scale must be non-negative")

    error_scale = scale * (dx ** accuracy)
    reduction_if_halved = 2 ** accuracy
    return {
        "error_scale": error_scale,
        "order": accuracy,
        "reduction_if_halved": reduction_if_halved,
    }


def parse_args() -> argparse.Namespace:
    parser = argparse.ArgumentParser(
        description="Estimate truncation error scaling for a scheme.",
        formatter_class=argparse.ArgumentDefaultsHelpFormatter,
    )
    parser.add_argument("--dx", type=float, required=True, help="Grid spacing")
    parser.add_argument("--accuracy", type=int, required=True, help="Scheme order")
    parser.add_argument(
        "--scale",
        type=float,
        default=1.0,
        help="Scale of higher derivative (C term)",
    )
    parser.add_argument("--json", action="store_true", help="Emit JSON output")
    return parser.parse_args()


def main() -> None:
    args = parse_args()
    try:
        result = estimate_truncation_error(args.dx, args.accuracy, args.scale)
    except ValueError as exc:
        print(str(exc), file=sys.stderr)
        sys.exit(2)

    payload = {
        "inputs": {
            "dx": args.dx,
            "accuracy": args.accuracy,
            "scale": args.scale,
        },
        "results": result,
    }

    if args.json:
        print(json.dumps(payload, indent=2, sort_keys=True))
        return

    print("Truncation error estimate")
    print(f"  error_scale: {result['error_scale']:.6g}")
    print(f"  order: {result['order']}")
    print(f"  reduction_if_halved: {result['reduction_if_halved']}")


if __name__ == "__main__":
    main()

```

### references/stencil_catalog.md

```markdown
# Stencil Catalog

Comprehensive reference for finite difference stencil coefficients.

## First Derivative Stencils

### Central Differences

Symmetric stencils, no directional bias.

**2nd Order (3-point):**
```
Offsets: [-1, 0, +1]
Coefficients: [-1/2, 0, 1/2] / dx

f'(x) ≈ (f(x+dx) - f(x-dx)) / (2dx)
Error: O(dx²)
```

**4th Order (5-point):**
```
Offsets: [-2, -1, 0, +1, +2]
Coefficients: [1/12, -2/3, 0, 2/3, -1/12] / dx

f'(x) ≈ (-f(x+2dx) + 8f(x+dx) - 8f(x-dx) + f(x-2dx)) / (12dx)
Error: O(dx⁴)
```

**6th Order (7-point):**
```
Offsets: [-3, -2, -1, 0, +1, +2, +3]
Coefficients: [-1/60, 3/20, -3/4, 0, 3/4, -3/20, 1/60] / dx
Error: O(dx⁶)
```

### Forward Differences (One-Sided)

For left boundaries or upwind schemes.

**1st Order (2-point):**
```
Offsets: [0, +1]
Coefficients: [-1, 1] / dx

f'(x) ≈ (f(x+dx) - f(x)) / dx
Error: O(dx)
```

**2nd Order (3-point):**
```
Offsets: [0, +1, +2]
Coefficients: [-3/2, 2, -1/2] / dx

f'(x) ≈ (-3f(x) + 4f(x+dx) - f(x+2dx)) / (2dx)
Error: O(dx²)
```

**3rd Order (4-point):**
```
Offsets: [0, +1, +2, +3]
Coefficients: [-11/6, 3, -3/2, 1/3] / dx
Error: O(dx³)
```

**4th Order (5-point):**
```
Offsets: [0, +1, +2, +3, +4]
Coefficients: [-25/12, 4, -3, 4/3, -1/4] / dx
Error: O(dx⁴)
```

### Backward Differences (One-Sided)

For right boundaries (mirror of forward).

**1st Order (2-point):**
```
Offsets: [-1, 0]
Coefficients: [-1, 1] / dx
```

**2nd Order (3-point):**
```
Offsets: [-2, -1, 0]
Coefficients: [1/2, -2, 3/2] / dx
```

## Second Derivative Stencils

### Central Differences

**2nd Order (3-point):**
```
Offsets: [-1, 0, +1]
Coefficients: [1, -2, 1] / dx²

f''(x) ≈ (f(x+dx) - 2f(x) + f(x-dx)) / dx²
Error: O(dx²)
```

**4th Order (5-point):**
```
Offsets: [-2, -1, 0, +1, +2]
Coefficients: [-1/12, 4/3, -5/2, 4/3, -1/12] / dx²

f''(x) ≈ (-f(x+2dx) + 16f(x+dx) - 30f(x) + 16f(x-dx) - f(x-2dx)) / (12dx²)
Error: O(dx⁴)
```

**6th Order (7-point):**
```
Offsets: [-3, -2, -1, 0, +1, +2, +3]
Coefficients: [1/90, -3/20, 3/2, -49/18, 3/2, -3/20, 1/90] / dx²
Error: O(dx⁶)
```

### One-Sided Second Derivative

**2nd Order Forward (4-point):**
```
Offsets: [0, +1, +2, +3]
Coefficients: [2, -5, 4, -1] / dx²
Error: O(dx²)
```

**2nd Order Backward (4-point):**
```
Offsets: [-3, -2, -1, 0]
Coefficients: [-1, 4, -5, 2] / dx²
Error: O(dx²)
```

## Higher Derivative Stencils

### Third Derivative

**2nd Order Central (5-point):**
```
Offsets: [-2, -1, 0, +1, +2]
Coefficients: [-1/2, 1, 0, -1, 1/2] / dx³
Error: O(dx²)
```

**4th Order Central (7-point):**
```
Offsets: [-3, -2, -1, 0, +1, +2, +3]
Coefficients: [1/8, -1, 13/8, 0, -13/8, 1, -1/8] / dx³
Error: O(dx⁴)
```

### Fourth Derivative

**2nd Order Central (5-point):**
```
Offsets: [-2, -1, 0, +1, +2]
Coefficients: [1, -4, 6, -4, 1] / dx⁴
Error: O(dx²)
```

**4th Order Central (7-point):**
```
Offsets: [-3, -2, -1, 0, +1, +2, +3]
Coefficients: [-1/6, 2, -13/2, 28/3, -13/2, 2, -1/6] / dx⁴
Error: O(dx⁴)
```

## Mixed and Cross Derivatives

### First Cross Derivative (∂²f/∂x∂y)

**2nd Order (4-point):**
```
f_xy ≈ (f(x+dx,y+dy) - f(x+dx,y-dy) - f(x-dx,y+dy) + f(x-dx,y-dy)) / (4 dx dy)

Stencil (in 2D):
    [+1]      [-1]
       [0]
    [-1]      [+1]
```

**4th Order:**
```
Uses 12 points in a wider pattern
Error: O(dx⁴ + dy⁴)
```

## Compact (Padé) Schemes

Higher accuracy with smaller stencils by solving implicit system.

### 4th Order Compact First Derivative

```
α f'_{i-1} + f'_i + α f'_{i+1} = (a/dx)(f_{i+1} - f_{i-1})

where α = 1/4, a = 3/4
```

Requires solving tridiagonal system for f'.

### 6th Order Compact First Derivative

```
α f'_{i-1} + f'_i + α f'_{i+1} = (a/dx)(f_{i+1} - f_{i-1}) + (b/dx)(f_{i+2} - f_{i-2})

where α = 1/3, a = 14/9, b = 1/9
```

### Advantages of Compact Schemes

| Property | Explicit | Compact |
|----------|----------|---------|
| Stencil width | Wide | Narrow |
| Operations | O(n) | O(n) but solve required |
| Spectral resolution | Good | Excellent |
| Boundary handling | Simple | Requires closure |

## Upwind Schemes

For advection-dominated problems with u > 0.

### First-Order Upwind

```
u > 0: f'(x) ≈ (f(x) - f(x-dx)) / dx  (backward)
u < 0: f'(x) ≈ (f(x+dx) - f(x)) / dx  (forward)
```

Stable but highly diffusive.

### Second-Order Upwind

```
u > 0: f'(x) ≈ (3f(x) - 4f(x-dx) + f(x-2dx)) / (2dx)
u < 0: f'(x) ≈ (-3f(x) + 4f(x+dx) - f(x+2dx)) / (2dx)
```

Less diffusive, may oscillate.

### Third-Order Upwind (QUICK)

```
f'(x) ≈ (2f(x+dx) + 3f(x) - 6f(x-dx) + f(x-2dx)) / (6dx)
```

Good balance of accuracy and stability.

## Stencil Properties Summary

### First Derivative

| Order | Points | Width | Coefficients Sum |
|-------|--------|-------|------------------|
| 2 central | 3 | 2dx | 0 |
| 4 central | 5 | 4dx | 0 |
| 6 central | 7 | 6dx | 0 |
| 1 one-sided | 2 | dx | 0 |
| 2 one-sided | 3 | 2dx | 0 |

### Second Derivative

| Order | Points | Width | Coefficients Sum |
|-------|--------|-------|------------------|
| 2 central | 3 | 2dx | 0 |
| 4 central | 5 | 4dx | 0 |
| 6 central | 7 | 6dx | 0 |

### Stability Properties

| Scheme | Dispersion | Dissipation |
|--------|------------|-------------|
| Central | Low | None |
| Upwind | Moderate | High |
| Compact | Very low | None |
| WENO | Low | Adaptive |

## Implementation Notes

### Applying Stencils

```python
def apply_stencil(f, coeffs, offsets, dx):
    """Apply finite difference stencil."""
    result = np.zeros_like(f)
    for c, o in zip(coeffs, offsets):
        result += c * np.roll(f, -o)
    return result / dx

# Example: 4th-order central first derivative
coeffs = [1/12, -2/3, 0, 2/3, -1/12]
offsets = [-2, -1, 0, 1, 2]
df_dx = apply_stencil(f, coeffs, offsets, dx)
```

### Boundary Treatment

Interior stencils don't work at boundaries:

```
For central 4th-order (needs ±2 points):
  x = 0: use forward stencil
  x = dx: use biased stencil
  x = 2dx to x = (n-3)dx: central stencil
  x = (n-2)dx: use biased stencil
  x = (n-1)dx: use backward stencil
```

### Verification

Taylor expansion verification:

```
f(x+dx) = f(x) + dx f'(x) + dx²/2 f''(x) + dx³/6 f'''(x) + ...

Substitute into stencil, verify leading error term.
```

Grid refinement verification:

```
Run with dx, dx/2, dx/4
Error should decrease as dx^p for order p
```

```

### references/boundary_handling.md

```markdown
# Boundary Handling

Comprehensive guide for implementing finite difference stencils at domain boundaries.

## The Boundary Problem

Interior stencils require points outside the domain:

```
Central 5-point stencil at x = dx:
  Needs: x - 2dx, x - dx, x, x + dx, x + 2dx
  But: x - 2dx = -dx is outside domain!
```

## One-Sided Stencils

### Reducing Stencil Order

Use lower-order one-sided stencils near boundaries:

```
Left boundary (x = 0):
  Point 0: 1st-order forward
  Point 1: 2nd-order forward (if needed)
  Points 2+: Central stencil

Right boundary (x = L):
  Point n-1: 1st-order backward
  Point n-2: 2nd-order backward (if needed)
  Points <n-2: Central stencil
```

### Matching Interior Order

To maintain overall accuracy, use one-sided stencils of same order as interior:

```
Interior: 4th-order central (5 points)
Boundary: 4th-order one-sided (5 points)

First derivative at x = 0:
  f'(0) ≈ (-25f_0 + 48f_1 - 36f_2 + 16f_3 - 3f_4) / (12dx)
```

### Biased Stencils

Use asymmetric stencils that include boundary:

```
At x = dx (one point from left boundary):
  Offsets: [-1, 0, +1, +2, +3]
  A 4th-order stencil shifted to fit

Coefficients (for f'):
  [-1/4, -5/6, 3/2, -1/2, 1/12] / dx
```

## Ghost Cell Method

### Concept

Add fictitious points outside domain:

```
Physical domain: x = 0, dx, 2dx, ..., (n-1)dx
Ghost cells: x = -dx, -2dx (left), x = n*dx, (n+1)*dx (right)

Now central stencils work at all interior points!
```

### Filling Ghost Cells

The ghost values must be computed from boundary conditions.

**Dirichlet BC:** f(0) = g

```
For 2nd-order central derivative at x = dx:
  Need f(-dx)
  Use symmetry: f(-dx) = 2g - f(dx) (linear extrapolation)
  Or: f(-dx) = -f(dx) + 2f(0) (reflect through boundary value)
```

**Neumann BC:** f'(0) = h

```
For 2nd-order:
  f'(0) ≈ (f(dx) - f(-dx)) / (2dx) = h
  → f(-dx) = f(dx) - 2dx*h

For 4th-order:
  Use higher-order approximation to get f(-dx), f(-2dx)
```

**Periodic BC:**

```
f(-dx) = f((n-1)*dx)
f(-2dx) = f((n-2)*dx)
f(n*dx) = f(0)
f((n+1)*dx) = f(dx)
```

### Ghost Cell Implementation

```python
def fill_ghosts(f, bc_type, bc_value, n_ghost):
    """Fill ghost cells based on boundary condition."""
    f_ext = np.zeros(len(f) + 2 * n_ghost)
    f_ext[n_ghost:-n_ghost] = f

    if bc_type == 'dirichlet':
        # Left: reflect through boundary value
        for i in range(n_ghost):
            f_ext[n_ghost - 1 - i] = 2 * bc_value[0] - f[i]
        # Right: reflect through boundary value
        for i in range(n_ghost):
            f_ext[-n_ghost + i] = 2 * bc_value[1] - f[-1 - i]

    elif bc_type == 'neumann':
        # Left: extrapolate with derivative
        for i in range(n_ghost):
            f_ext[n_ghost - 1 - i] = f[i] - 2 * (i + 1) * dx * bc_value[0]
        # Right: extrapolate with derivative
        for i in range(n_ghost):
            f_ext[-n_ghost + i] = f[-1 - i] + 2 * (i + 1) * dx * bc_value[1]

    elif bc_type == 'periodic':
        f_ext[:n_ghost] = f[-n_ghost:]
        f_ext[-n_ghost:] = f[:n_ghost]

    return f_ext
```

## Boundary Condition Types

### Dirichlet (Specified Value)

```
f(boundary) = g(t)

Direct substitution:
  f_0 = g  (known, not part of solve)

For Laplacian discretization:
  (f_1 - 2f_0 + f_{-1}) / dx² at x = dx
  becomes: (f_1 - 2g + f_{-1}) / dx²
  where f_{-1} is ghost (computed for consistency)
```

### Neumann (Specified Derivative)

```
f'(boundary) = h(t)

Using one-sided difference at x = 0:
  (f_1 - f_0) / dx = h → f_0 = f_1 - dx * h

Or using ghost cell:
  (f_1 - f_{-1}) / (2dx) = h → f_{-1} = f_1 - 2dx * h
```

### Robin (Mixed)

```
a * f(boundary) + b * f'(boundary) = c

Combine Dirichlet and Neumann treatment.
```

### Periodic

```
f(0) = f(L), f'(0) = f'(L), etc.

Simply wrap indices:
  f_{-1} = f_{n-1}
  f_{n} = f_0
```

### Symmetry

```
Even symmetry at x = 0: f(-x) = f(x)
  → f_{-1} = f_1, f_{-2} = f_2

Odd symmetry at x = 0: f(-x) = -f(x)
  → f_{-1} = -f_1, f_{-2} = -f_2
```

## Accuracy Considerations

### Order Reduction at Boundaries

Using lower-order BC treatment reduces global accuracy:

| Interior Order | BC Treatment | Global Order |
|----------------|--------------|--------------|
| 4 | 1st-order | 1 |
| 4 | 2nd-order | 2 |
| 4 | 3rd-order | 3 |
| 4 | 4th-order | 4 |

**Rule of thumb:** BC treatment should be at least (interior order - 1).

### Consistency Requirements

For 4th-order Laplacian with Neumann BC:

```
Need: f'(0) to O(dx³) at minimum
      f(-dx), f(-2dx) consistent with this

Standard approach:
  1. Use 4th-order formula for f'(0) = h
  2. Solve for ghost values
  3. Apply interior stencil
```

## Special Cases

### Corner Boundaries (2D/3D)

At corners, multiple boundaries meet:

```
Corner at (0, 0):
  Need ghost values at (-dx, 0), (0, -dy), (-dx, -dy)

Strategy 1: Edge then corner
  Fill edge ghosts first, then extrapolate corners

Strategy 2: Direct corner treatment
  Use corner-specific stencils
```

### Moving Boundaries

For domains with moving interfaces:

```
1. Identify boundary location (between grid points)
2. Use immersed boundary or ghost fluid methods
3. Extrapolate/interpolate to find ghost values
```

### Non-Uniform Grids

Ghost cell spacing may differ from interior:

```
If dx varies:
  Use variable-coefficient stencils near boundary
  Or: Map to uniform computational coordinates
```

## Implementation Patterns

### Modular BC Application

```python
class BoundaryCondition:
    def apply(self, f, dx):
        raise NotImplementedError

class DirichletBC(BoundaryCondition):
    def __init__(self, value):
        self.value = value

    def apply(self, f, dx):
        f[0] = self.value  # Direct enforcement

class NeumannBC(BoundaryCondition):
    def __init__(self, gradient):
        self.gradient = gradient

    def apply(self, f, dx):
        # Second-order one-sided
        f[0] = f[1] - dx * self.gradient
```

### Matrix Form for Implicit

When forming Ax = b:

```
Dirichlet at i = 0:
  A[0, 0] = 1
  A[0, 1:] = 0
  b[0] = g

Neumann at i = 0 (2nd order):
  A[0, 0] = -1/dx
  A[0, 1] = 1/dx
  b[0] = h
```

## Verification

### Manufactured Solutions

Test with known solution that satisfies BCs:

```
f(x) = sin(πx) on [0, 1]
Dirichlet: f(0) = 0, f(1) = 0 ✓

f(x) = cos(πx) on [0, 1]
Neumann: f'(0) = 0, f'(1) = 0 ✓
```

### Convergence Test

Refine grid, verify expected order:

```
dx, dx/2, dx/4:
  Error should decrease as dx^p

If not: BC treatment is limiting accuracy
```

### Symmetry Test

For symmetric problems:

```
If physics is symmetric about x = L/2:
  Solution should be symmetric
  Check: max|f(x) - f(L-x)| < tolerance
```

## Common Pitfalls

### Pitfall 1: Inconsistent Ghost Values

```
Problem: Ghost values computed at wrong order
Result: Loss of accuracy, may look like noise near boundary
Fix: Ensure ghost fill matches interior stencil order
```

### Pitfall 2: Wrong Boundary Location

```
Problem: BC applied at wrong grid point (cell center vs face)
Result: O(dx) error in solution
Fix: Be clear about grid convention
```

### Pitfall 3: Time-Dependent BC Lag

```
Problem: BC evaluated at old time in implicit scheme
Result: Loss of time accuracy
Fix: Use BC at new time (may require iteration)
```

### Pitfall 4: Missing Corner Treatment

```
Problem: Corners filled by extrapolation that doesn't match physics
Result: Artifacts at corners
Fix: Proper corner BC (often requires special handling)
```

```

### references/scheme_selection.md

```markdown
# Scheme Selection

Comprehensive guide for choosing numerical differentiation schemes.

## Scheme Categories

### Finite Difference (FD)

Approximate derivatives using point values on a grid.

**Characteristics:**
- Simple implementation
- Works on structured grids
- Accuracy depends on stencil width

**Best for:**
- Smooth solutions
- Structured rectangular domains
- Moderate accuracy requirements

### Finite Volume (FV)

Discretize integral conservation laws.

**Characteristics:**
- Inherently conservative
- Natural for conservation laws
- Works on unstructured meshes

**Best for:**
- Discontinuous solutions
- Conservation-critical problems
- Complex geometries

### Finite Element (FE)

Approximate solution as sum of basis functions.

**Characteristics:**
- Flexible geometry handling
- Natural for variational problems
- Higher-order possible

**Best for:**
- Complex geometries
- Multi-physics problems
- Adaptive refinement

### Spectral Methods

Represent solution in global basis (Fourier, Chebyshev).

**Characteristics:**
- Exponential accuracy for smooth problems
- Global operations (FFT)
- Sensitive to discontinuities

**Best for:**
- Periodic domains
- Very smooth solutions
- High accuracy requirements

## Decision Flowchart

```
What is the solution character?
│
├── Smooth everywhere
│   ├── Periodic domain? → Spectral (FFT)
│   ├── Simple geometry? → High-order FD
│   └── Complex geometry? → FE or FD on curvilinear
│
├── Contains discontinuities (shocks, interfaces)
│   ├── Conservation critical? → FV with limiters
│   └── Conservation not critical? → FD with WENO
│
└── Mixed (smooth regions + discontinuities)
    ├── Localized discontinuity? → Hybrid FD + shock capturing
    └── Many discontinuities? → FV with AMR
```

## Detailed Comparisons

### FD vs FV

| Aspect | Finite Difference | Finite Volume |
|--------|-------------------|---------------|
| Formulation | Pointwise derivatives | Cell-averaged fluxes |
| Conservation | Not automatic | Built-in |
| Shocks | Needs limiting | Natural handling |
| Grids | Structured preferred | Any mesh |
| Implementation | Simpler | More complex |
| Accuracy order | Easy high-order | High-order harder |

### Central vs Upwind FD

| Aspect | Central | Upwind |
|--------|---------|--------|
| Bias | None | Flow direction |
| Accuracy | Higher order | Lower order |
| Dispersion | Low | Higher |
| Dissipation | None | High |
| Stability | Needs care | Naturally stable |
| Use for | Diffusion | Advection |

### Explicit vs Implicit

| Aspect | Explicit FD | Implicit FD |
|--------|-------------|-------------|
| Time step | Limited by CFL | Not limited |
| Cost per step | Low | High (solve) |
| Stiffness | Cannot handle | Can handle |
| Parallelism | Easy | Harder |
| Implementation | Simple | Complex |

## Specific Applications

### Diffusion Equation

```
∂u/∂t = D∇²u
```

**Recommended:**
- Spatial: Central FD (2nd or 4th order)
- Time: Implicit (Crank-Nicolson) or explicit RK

**Avoid:**
- Upwind (diffusion has no direction)
- First-order time (too much numerical diffusion)

### Advection Equation

```
∂u/∂t + v·∇u = 0
```

**Recommended:**
- Upwind FD (stable, diffusive)
- Central + artificial viscosity
- WENO for sharp features
- Semi-Lagrangian for large CFL

**Avoid:**
- Pure central difference (unstable without stabilization)

### Wave Equation

```
∂²u/∂t² = c²∇²u
```

**Recommended:**
- Central FD in space
- Leap-frog or RK in time
- Symplectic for long-time accuracy

**Avoid:**
- Upwind (too diffusive)
- Low-order time stepping (phase errors)

### Convection-Diffusion

```
∂u/∂t + v·∇u = D∇²u
```

**Recommended:**
- Peclet-dependent blending
- SUPG (Streamline Upwind Petrov-Galerkin)
- High Pe: upwind-biased
- Low Pe: central

### Nonlinear Conservation Laws

```
∂u/∂t + ∇·f(u) = 0
```

**Recommended:**
- FV with Riemann solvers
- ENO/WENO
- Flux limiters (minmod, Van Leer, etc.)

### Phase-Field Equations

```
∂φ/∂t = -M δF/δφ
```

**Recommended:**
- Central FD for Laplacian
- Conservative form for Cahn-Hilliard
- High-order for interface resolution

## Order Selection

### When to Use Low Order (2nd)

- Initial development/debugging
- Solutions are smooth but not analytic
- Memory is limiting factor
- Complex boundary conditions
- Unstructured or adaptive meshes

### When to Use High Order (4th+)

- Very smooth solutions
- High accuracy needed (validation)
- Uniform structured grids
- Long-time integration (error accumulation)
- Wave propagation (reduce dispersion)

### When to Use Spectral

- Periodic domains
- Smooth, analytic solutions
- Turbulence (DNS)
- Maximum accuracy for given resolution

## Stability-Accuracy Trade-offs

### Accuracy vs Stability

| Choice | Accuracy | Stability |
|--------|----------|-----------|
| Central + fine dt | High | OK with small dt |
| Upwind | Low | Very stable |
| Central + stabilization | Medium | Stable |
| WENO | High | Stable |
| Implicit central | High | Unconditional |

### Computational Cost

| Scheme | Operations/point | Memory |
|--------|------------------|--------|
| 2nd-order central | 3 | O(1) |
| 4th-order central | 5 | O(1) |
| Compact 6th-order | 3 + solve | O(n) |
| Spectral | O(n log n) FFT | O(n) |
| WENO5 | ~15 | O(1) |

## Grid Considerations

### Uniform Grids

All schemes work well. Use highest practical order.

### Stretched Grids

Transform to computational coordinates:

```
x → ξ (uniform)
dx/dξ varies

∂f/∂x = (∂f/∂ξ) / (∂x/∂ξ)
```

Or use variable-coefficient stencils.

### Unstructured Grids

Limited options:
- FV (natural)
- FE (natural)
- FD on stencils (least-squares fitting)

### Adaptive Meshes

Need:
- Conservative interpolation at refinement boundaries
- Consistent stencils across levels
- Flux matching for conservation

## Scheme Selection Checklist

- [ ] Identified solution character (smooth/discontinuous)
- [ ] Determined conservation requirements
- [ ] Assessed grid structure (structured/unstructured)
- [ ] Balanced accuracy vs cost
- [ ] Considered stability constraints
- [ ] Checked boundary condition compatibility
- [ ] Verified with manufactured solution
- [ ] Performed grid convergence study

## Quick Reference Table

| Problem | Smooth | Periodic | Conservation | Recommended |
|---------|--------|----------|--------------|-------------|
| Diffusion | Yes | No | No | Central FD |
| Diffusion | Yes | Yes | No | Spectral |
| Advection | Yes | Any | No | Central + stab |
| Advection | No | Any | Yes | FV + limiter |
| Wave | Yes | No | No | Central + symplectic |
| Wave | Yes | Yes | No | Spectral |
| Navier-Stokes | Mixed | No | Yes | FV or FD + projection |
| Phase-field | Yes | No | Maybe | High-order FD |

```

### references/error_guidance.md

```markdown
# Truncation Error Guidance

Comprehensive guide for understanding and controlling discretization errors.

## Truncation Error Fundamentals

### Definition

Truncation error is the difference between the continuous operator and its discrete approximation:

```
τ = L[u] - L_h[u_h]

where:
  L[u] = continuous differential operator applied to exact solution
  L_h[u_h] = discrete operator applied to grid values
```

### Order of Accuracy

A scheme is pth-order accurate if:

```
τ = O(dx^p)  as  dx → 0

Meaning: τ ≤ C × dx^p for some constant C
```

| Order p | Error Reduction with dx/2 |
|---------|---------------------------|
| 1 | 2× |
| 2 | 4× |
| 3 | 8× |
| 4 | 16× |
| 6 | 64× |

## Error Estimation

### Taylor Expansion Method

For 2nd-order central first derivative:

```
f'(x) ≈ (f(x+dx) - f(x-dx)) / (2dx)

Taylor expand:
f(x+dx) = f(x) + dx f'(x) + (dx²/2) f''(x) + (dx³/6) f'''(x) + ...
f(x-dx) = f(x) - dx f'(x) + (dx²/2) f''(x) - (dx³/6) f'''(x) + ...

Subtract: f(x+dx) - f(x-dx) = 2dx f'(x) + (dx³/3) f'''(x) + ...

Divide: (f(x+dx) - f(x-dx)) / (2dx) = f'(x) + (dx²/6) f'''(x) + ...

Error = (dx²/6) f'''(x) + O(dx⁴)
```

### Leading Error Term

The coefficient of the leading error term matters:

| Scheme | Order | Leading Coefficient |
|--------|-------|---------------------|
| Central d/dx | 2 | dx²/6 × f''' |
| Central d²/dx² | 2 | dx²/12 × f'''' |
| 4th-order d/dx | 4 | dx⁴/30 × f''''' |
| Upwind d/dx | 1 | dx/2 × f'' |

## Practical Error Estimation

### Richardson Extrapolation

Compare solutions at different resolutions:

```
Solution with dx: u_h
Solution with dx/2: u_{h/2}

For pth-order scheme:
  u_h = u_exact + C × dx^p + O(dx^{p+1})
  u_{h/2} = u_exact + C × (dx/2)^p + O(dx^{p+1})

Subtract:
  u_h - u_{h/2} ≈ C × dx^p × (1 - 1/2^p)

Error estimate:
  |u_h - u_exact| ≈ |u_h - u_{h/2}| / (2^p - 1)
```

### Order Verification

```
e_h = ||u_h - u_exact||
e_{h/2} = ||u_{h/2} - u_exact||

Observed order: p_obs = log(e_h / e_{h/2}) / log(2)

Should match theoretical order p.
```

### Convergence Study

```python
def convergence_study(exact, solve, dx_list):
    """Compute errors and observed order."""
    errors = []
    for dx in dx_list:
        u = solve(dx)
        u_exact = exact(grid(dx))
        errors.append(np.linalg.norm(u - u_exact))

    orders = []
    for i in range(len(errors) - 1):
        ratio = errors[i] / errors[i+1]
        dx_ratio = dx_list[i] / dx_list[i+1]
        orders.append(np.log(ratio) / np.log(dx_ratio))

    return errors, orders
```

## Resolution Requirements

### Points Per Feature

Minimum grid points to resolve a feature:

| Feature | Min Points | Recommended |
|---------|------------|-------------|
| Wavelength | 5 | 10-20 |
| Boundary layer | 5 | 10 |
| Phase-field interface | 3 | 5-10 |
| Shock width | 2-3 | 5+ |
| Vortex | 5 | 10-20 |

### Nyquist Limit

For oscillatory solutions:

```
Minimum: dx < λ/2 (Nyquist)
Practical: dx < λ/10 to λ/20 for accuracy
```

### Error vs Resolution

```
Error ≈ C × (dx/L_feature)^p

where L_feature is the smallest feature size.

Example: 2nd-order scheme, L = 0.1, dx = 0.01
  Error ≈ C × (0.01/0.1)² = C × 0.01 = 1% of C

To get 0.01%: dx = 0.001
```

## Common Error Sources

### Spatial Discretization Error

```
Source: Finite difference approximation
Scales as: O(dx^p)
Control: Refine grid or use higher-order scheme
```

### Temporal Discretization Error

```
Source: Time integration scheme
Scales as: O(dt^q) for qth-order time scheme
Control: Reduce dt or use higher-order integrator
```

### Roundoff Error

```
Source: Floating-point arithmetic
Scales as: O(ε_machine / dx^p) for pth derivative
Warning: Can dominate for very fine grids!
```

### Iteration Error

```
Source: Incomplete convergence of iterative solvers
Control: Tighten solver tolerance
Rule: Solver tolerance < discretization error
```

## Error Balancing

### Time-Space Error Balance

For optimal efficiency:

```
Temporal error ≈ Spatial error

dt^q ≈ dx^p

For explicit parabolic (CFL): dt ~ dx²
  If time is O(dt), space is O(dx²) → balanced for 2nd-order

For explicit hyperbolic (CFL): dt ~ dx
  If time is O(dt), space is O(dx²) → space limits accuracy
```

### Cost-Accuracy Trade-off

| Scheme | Error | Cost | Efficiency |
|--------|-------|------|------------|
| 2nd-order | O(dx²) | O(1/dx^d) | Baseline |
| 4th-order | O(dx⁴) | O(1/dx^d) | Better for smooth |
| 6th-order | O(dx⁶) | O(1/dx^d) | Best for very smooth |

For same accuracy, higher order needs fewer points.

## Validation Strategies

### Manufactured Solutions

1. Choose smooth function: u(x,t) = sin(πx)cos(t)
2. Substitute into PDE: get source term
3. Solve with source
4. Compare to known u

**Advantages:**
- Known exact solution
- Tests all error sources
- Verifies implementation

### Grid Convergence

1. Solve on dx, dx/2, dx/4
2. Compute errors vs fine grid or analytical
3. Verify expected order

```
If order drops:
  - Bug in implementation
  - Solution not smooth enough
  - Boundary conditions limiting
  - Grid not fine enough yet
```

### Conservation Check

For conservation laws:

```
∫ u dx at t=0 should equal ∫ u dx at t=T (up to BCs)

If not conserved:
  - Non-conservative scheme
  - Discretization error
  - Boundary flux errors
```

## Troubleshooting

### Order Lower Than Expected

| Symptom | Possible Cause | Fix |
|---------|----------------|-----|
| p_obs = 1 | Boundary limiting | Higher-order BC |
| p_obs = 0 | Solution not smooth | Check physics |
| p_obs varies | Pre-asymptotic | Finer grids |
| p_obs > expected | Super-convergence | Lucky, verify |

### Error Not Decreasing

| Symptom | Possible Cause | Fix |
|---------|----------------|-----|
| Error constant | Roundoff limit | Coarsen grid |
| Error increases | Instability | Check dt, scheme |
| Error oscillates | Non-monotone | Different norms |

### Unexpected Error Growth in Time

| Symptom | Possible Cause | Fix |
|---------|----------------|-----|
| Linear growth | Dispersion error | Higher order |
| Exponential growth | Instability | Check CFL |
| Bounded but large | Phase error | Smaller dt |

## Error Norms

### Common Norms

| Norm | Formula | Emphasis |
|------|---------|----------|
| L∞ (max) | max\|e_i\| | Worst point |
| L2 (RMS) | sqrt(Σ e_i² / n) | Average |
| L1 (mean) | Σ \|e_i\| / n | Mild outliers |
| Weighted | Σ w_i e_i² | Specific regions |

### Norm Selection

- **L∞**: When worst case matters (stability)
- **L2**: General accuracy assessment
- **L1**: When outliers expected (discontinuities)
- **Weighted**: Focus on region of interest

## Quick Reference

### Error Scaling Summary

| dx | 2nd-order error | 4th-order error |
|----|-----------------|-----------------|
| 0.1 | 0.01 | 0.0001 |
| 0.05 | 0.0025 | 6.25e-6 |
| 0.01 | 0.0001 | 1e-8 |

### Target Accuracy Guidelines

| Application | Typical Accuracy | Grid Requirement |
|-------------|------------------|------------------|
| Engineering | 1% | ~10 points/feature |
| Research | 0.1% | ~20 points/feature |
| Validation | 0.01% | ~50 points/feature |
| Spectral accuracy | 1e-10 | Smooth + spectral method |

```