I need to populate the matrix (stored as an array of arrays) with some values. The matrix is a Jacobian for a simple diffusion problem and looks like this:
J(1,1) = 1, J(N,N)=0
and for 1<n<N:
J(n,n) = -2k/dx^2 - 2*c(n)
J(n,n-1)=J(n,n+1) = k/dx^2
the rest of the matrix entries are zeros.
So far I have this monstrosity:
(1 to: c size) collect: [ :n |
(1 to: c size) collect: [ :m |
n = 1 | (n = c size)
ifTrue: [ m = n ifTrue: [ 1.0 ] ifFalse: [ 0.0 ] ]
ifFalse: [ m = n
ifTrue: [ -2.0 * k / dx squared - (2.0 * (c at: n)) ]
ifFalse: [ m = (n-1) | (m = (n+1))
ifTrue: [ k / dx squared ]
ifFalse: [ 0.0 ] ] ]
] ]
Notice the nested "if-statements" (Smalltalk equivalents). This works. But, is there, perhaps, a more elegant way of doing the same thing? As it stands now, it is rather unreadable.
For readability's sake I would consider sacrificing the extra O(n) time and avoid IFs altogether (which just make it even faster...).
J(N,N) = 0
J(1,1) = 1
//and for 1<n<N:
J(n,n) = Y(n)
J(n,m-1) = J(n,m+1) = X
What this tells me is that the whole matrix looks something like this
( 1 X 0 0 0 )
( X Y X 0 0 )
( 0 X Y X 0 )
( 0 0 X Y X )
( 0 0 0 X 0 )
Which means that I can create the whole matrix with just zeros, and then change the diagonal and neighboring diagonals.
jNM := [ k / dx squared ].
jNN := [ :n | -2.0 * k / dx squared - (2.0 * (c at: n)) ].
n := c size.
m := Matrix
new: n
tabulate: [:i :j | 0 ].
(1 to: n - 1) do: [ :i |
m at: i at: i put: (jNN value: i).
m at: i + 1 at: i put: jnM value.
m at: i at: i + 1 put: jnM value.
].
m at: 1 at: 1 put: 1.
Note: I'm not familiar with the math behind this but the value for J(n,m-1) seems like a constant to me.
Note 2: I'm putting the values at i + 1 indexes, because I am starting at position 1;1, but you can start from the opposite direction and have i-1.