better documentation

README.md

@@ -1,3 +1,48 @@
 ## Option Price PDE
 
-This is a basic library for solving 1-dimensional PDEs of the form `f_t = alpha(x, t) f_x + sigma(x, t) f_xx + rho(x, t) f`.  It uses finite differences and implicit solutions using the Thomas algorithm.  It is tested on the "basic" Black Scholes equation (with constant coefficients via the transformation `x = log(S)`) and with a transformation that maps the infinite Black Scholes domain to a finite domain (via the transformation `x = S/(S+K)`)
+This is a basic library for solving 1-dimensional PDEs of the form `f_t = alpha(x, t) f_x + sigma(x, t) f_xx + rho(x, t) f`.  It uses finite differences and implicit solutions using the Thomas algorithm.  It is tested on the "basic" Black Scholes equation (with constant coefficients via the transformation `x = log(S)`) and with a transformation that maps the infinite Black Scholes domain to a finite domain (via the transformation `x = S/(S+K)`)
+
+
+## Example
+
+```rust
+extern crate black_scholes;
+extern crate option_price_pde;
+let strike = 4.5;
+let rate = 0.01;
+let sigma = 0.3;
+let maturity = 2.0;
+let t_lim: solve_second_order_pde::Boundary = solve_second_order_pde::Boundary {
+    min: 0.0,
+    max: maturity,
+};
+let x_lim = solve_second_order_pde::Boundary { min: 0.0, max: 1.0 };
+let size_t = 700;
+let size_x = 255;
+let result = solve_second_order_pde::solve_second_order_pde(
+    |_t, x| -rate * (1.0 - x),                                        // coefficient on f(x)
+    |_t, x| rate * (1.0 - x) * x,                                     // coefficient on f'(x)
+    |_t, x: f64| (1.0 - x) * (1.0 - x) * x * x * sigma * sigma * 0.5, // coefficient on f''(x)
+    |_t, x| (2.0 * x - 1.0),                                          // upper boundary
+    |_t, _x| 0.0,                                                     // lower boundary
+    |x| {                                                             // time boundary (t=0)
+        let v = 2.0 * x - 1.0;
+        if v > 0.0 {
+            v
+        } else {
+            0.0
+        }
+    },
+    700,
+    &t_lim,
+    255,
+    &x_lim,
+)
+.unwrap();
+let (_t, x, v) = result
+    .get_result_at_x_and_t(size_t - 1, ((size_x - 1) / 2) as usize)
+    .unwrap();
+let s = (x * strike) / (1.0 - x);
+let approx_bs = (s + strike) * v;
+assert!((approx_bs - black_scholes::call(s, strike, rate, sigma, maturity)).abs() < 0.001);
+```

src/lib.rs

@@ -80,7 +80,7 @@ pub fn solve_second_order_pde(
     fn_fn: impl Fn(f64, f64) -> f64, //for black scholes, this is a constant (-r)
     fn_first_deriv: impl Fn(f64, f64) -> f64, //for black scholes, this is r*S (or, in logs, a constant r-sigma*sigma/2)
     fn_second_deriv: impl Fn(f64, f64) -> f64, //for black scholes, this is sigma*sigma*S*S/2
-    boundary_up: impl Fn(f64, f64) -> f64,    //for black scholes, this is s-ke^{-rt}
+    boundary_up: impl Fn(f64, f64) -> f64,    //for black scholes, this is s-k
     boundary_down: impl Fn(f64, f64) -> f64,  //for black scholes, this is 0
     init_fn: impl Fn(f64) -> f64,             //for black scholes, this is (s-k)^+
     size_t: usize,
@@ -176,7 +176,7 @@ mod tests {
             |_t, _x| -rate,
             |_t, _x| rate - sigma * sigma * 0.5,
             |_t, _x: f64| sigma * sigma * 0.5,
-            |t, x| (stock * x.exp() - strike) * (-rate * t).exp(),
+            |_t, x| (stock * x.exp() - strike),
             |_t, _x| 0.0,
             |x| {
                 let v = stock * x.exp() - strike;
@@ -192,7 +192,7 @@ mod tests {
             &x_lim,
         )
         .unwrap();
-        assert_eq!(result.surface.len(), size_t); //starts at "0"
+        assert_eq!(result.surface.len(), size_t);
         let (_t, x, approx_bs) = result
             .get_result_at_x_and_t(size_t - 1, ((size_x - 1) / 2) as usize)
             .unwrap();
@@ -218,7 +218,7 @@ mod tests {
             |_t, x| -rate * (1.0 - x),
             |_t, x| rate * (1.0 - x) * x,
             |_t, x: f64| (1.0 - x) * (1.0 - x) * x * x * sigma * sigma * 0.5,
-            |_t, x| (2.0 * x - 1.0), // * (-rate * t).exp(),
+            |_t, x| (2.0 * x - 1.0),
             |_t, _x| 0.0,
             |x| {
                 let v = 2.0 * x - 1.0;