<jupyter_start><jupyter_text>Introduction This is a short example of how to use the [SUNDIALS](https://computing.llnl.gov/projects/sundials) to solvea simple partial differential equation in Haskell via the [hmatrix-sundials](http://hackage.haskell.org/package/hmatrix-sundials)library. The example is taken from the C examples that come with the SUNDIALS source.  A Bit of Theory Consider the heat equation. With one spatial dimension we have:$$\frac{\partial u}{\partial t}=k \frac{\partial^{2} u}{\partial x^{2}}+f$$where $k$ is some diffusion coefficient. Let's take $t \in [0,10]$ and $x \in [0,1]$ with initial condition $u(0, x)=0$. Stationary boundary conditions$$\frac{\partial u}{\partial t}(t, 0)=\frac{\partial u}{\partial t}(t, 1)=0$$ And a point heat source$$f(t, x)=\left\{\begin{array}{ll}{1} & {\text { if } x=1 / 2} \\ {0} & {\text { otherwise }}\end{array}\right.$$ and we can discretize over this spatial dimension using:$$u_{x x}=\frac{u_{j+1}-2 u_{j}+u_{j-1}}{\Delta x^{2}}$$where$$u_{j}(t) \triangleq u\left(t, x_{j}\right), \quad x_{j} \triangleq j \Delta x, \quad 0 \leq j \leq n+1$$ We can write this in matrix form$$\dot{u}_i = \sum_0^{n+1} A_{i\,j} u_j + B_i, \quad 0 \leq i \leq n+1$$where$$\begin{aligned}A_{0\,j}     = 0, & \quad 0 \leq j \leq n+1, & \text{boundary condition} \\A_{i\,i-1}   = 1  &                          &                           \\A_{i\,i}     = 2  &                          &                           \\A_{i\,i+1}   = 1  &                          &                           \\A_{{n+1}\,j} = 0, & \quad 0 \leq j \leq n+1, & \text{boundary condition} \\A_{i\,j}     = 0  & \quad \text{otherwise}   &                           \\\end{aligned}$$ Converting this to a system of ODEs is straightforward, for example$$\begin{bmatrix}\dot{u_0} \\\dot{u_1} \\\dot{u_2} \\\dot{u_3} \\\dot{u_4}\end{bmatrix}=\begin{bmatrix}0 &  0 &  0 &  0 & 0 \\1 & -2 &  1 &  0 & 0 \\0 &  1 & -2 &  1 & 0 \\0 &  0 &  1 & -2 & 1 \\0 &  0 &  0 &  0 & 0\end{bmatrix}\begin{bmatrix}u_0 \\u_1 \\u_2 \\u_3 \\u_4\end{bmatrix}+\begin{bmatrix}f_0 \\f_1 \\f_2 \\f_3 \\f_4\end{bmatrix}$$where $f_j \triangleq f(t, x_j)$.  Implementation<jupyter_code>import           Numeric.Sundials.ARKode.ODE
import           Numeric.LinearAlgebra<jupyter_output><empty_output><jupyter_text>Spatial mesh size<jupyter_code>bigN :: Int
bigN = 201<jupyter_output><empty_output><jupyter_text>Heat conductivity<jupyter_code>k :: Double
k = 0.5<jupyter_output><empty_output><jupyter_text>Mesh spacing<jupyter_code>deltaX :: Double
deltaX = 1.0 / (fromIntegral bigN - 1)
c1, c2 :: Double
c1 = k / deltaX / deltaX
c2 = (-2.0) * k / deltaX / deltaX<jupyter_output><empty_output><jupyter_text>Initial time<jupyter_code>t0 :: Double
t0 = 0.0<jupyter_output><empty_output><jupyter_text>Final time<jupyter_code>tf :: Double
tf =1.0<jupyter_output><empty_output><jupyter_text>Total number of output times<jupyter_code>bigNt :: Int
bigNt = 10<jupyter_output><empty_output><jupyter_text>Timesteps<jupyter_code>deltaT :: Double
deltaT = (tf - t0) / (fromIntegral bigNt)<jupyter_output><empty_output><jupyter_text>Relative tolerance<jupyter_code>rtol :: Double
rtol = 1.0e-6<jupyter_output><empty_output><jupyter_text>Absolute tolerance<jupyter_code>atol :: Double
atol = 1.0e-10
bigA :: Matrix Double
bigA = assoc (bigN, bigN) 0.0 [ ((i, j), f (i, j)) | i <- [0 .. bigN - 1]
                                                   , j <- [0 .. bigN - 1]
                        ]
 where
   f (i, j) | i     == 0        = 0.0    -- left boundary condition
            | i     == bigN - 1 = 0.0    -- right boundary condition
            | i     == j        = c2
            | i - 1 == j        = c1
            | i + 1 == j        = c1
            | otherwise         = 0.0
b :: Vector Double
b = assoc bigN 0.0 [ (iSource, 0.01 / deltaX) ]
  where
    iSource = bigN `div` 2<jupyter_output><empty_output><jupyter_text>Initial condition<jupyter_code>bigU0 :: Vector Double
bigU0 = assoc bigN 0.0 []<jupyter_output><empty_output><jupyter_text>Finally we get to solve the system of equations<jupyter_code>sol :: Matrix Double
sol = odeSolveV SDIRK_5_3_4' Nothing rtol atol (const bigU') bigU0 (vector (map (deltaT *) [0 .. 10]))
  where
    bigU' bigU = bigA #> bigU + b<jupyter_output><empty_output><jupyter_text>We can draw the results on a chart using the [diagrams](http://hackage.haskell.org/package/diagrams)and [Chart](http://hackage.haskell.org/package/Chart) packages.<jupyter_code>import Diagrams.Prelude
import Graphics.Rendering.Chart
import Data.Default.Class
import Control.Lens
chart t l obss = toRenderable layout
  where

    actual x l c = plot_lines_values .~ [x]
                   $ plot_lines_style  . line_color .~ opaque c
                   $ plot_lines_title .~ l
                   $ plot_lines_style  . line_width .~ 1.0
                   $ def

    ls = map (\n -> "Time " ++ show (deltaT * fromIntegral n)) [0..]
    cs = cycle [blue, green, red]

    actuals' :: [PlotLines Double Double]
    actuals' = zipWith3 actual obss ls cs

    layout = layout_title .~ t
           $ layout_plots .~ map toPlot actuals'
           $ layout_y_axis . laxis_title .~ l
           $ layout_y_axis . laxis_override .~ axisGridHide
           $ layout_x_axis . laxis_title .~ "Time"
           $ layout_x_axis . laxis_override .~ axisGridHide
           $ def
chart "Simple 1D Heat Equation" "Temperature"
      (map (zip (map ((* deltaX) . fromIntegral) [0 .. bigN]) . toList) (toRows sol))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14976457805143065, 'fraction_numerical': 0.027526258601955813, 'mean_word_length': 2.657615894039735, 'pattern_counts': {'":': 0, '<': 55, '<?xml version=': 0, '>': 55, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '124770', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 2089, 'n_tokens_neox': 1972, 'n_words': 490, 'path': 'One Dimensional Heat Equation.ipynb', 'repo_name': 'idontgetoutmuch/Diffusions', 'size': '5522'}