我需要用一些值填充矩阵(存储为数组数组)。矩阵是一个简单的扩散问题的雅可比矩阵,看起来像这样:
J(1,1) = 1, J(N,N)=0
对于1<n<N
:
J(n,n) = -2k/dx^2 - 2*c(n)
J(n,n-1)=J(n,n+1) = k/dx^2
其余矩阵条目为零。
到目前为止,我有这个怪物:
(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 ] ] ]
] ]
注意嵌套的"if-statement"(Smalltalk等效项)。这行得通。但是,也许有一种更优雅的方式来做同样的事情?就目前而言,它是相当不可读的。
n := c size.
Matrix
new: n
tabulate: [:i :j | self jacobianAtRow: i column: j]
哪里
jacobianAtRow: i column: j
n := c size.
(i = 1 or: [i = n]) ifTrue: [^j = i ifTrue: [1.0] ifFalse [0.0]].
j = i ifTrue: [^-2.0 * k / dx squared - (2.0 * (c at: i))].
(j = (i - 1) or: [j = (i + 1)]) ifTrue: [^k / dx squared].
^0.0
基本上,一般的想法是这样的:每当你找到嵌套的ifs时,将这段代码分解成一个方法本身,并将嵌套转换为一个类似大小写的枚举,在每种可能性下都返回一个值。
为了可读性,我会考虑牺牲额外的O(n)
时间并完全避免 IF(这只会让它更快......
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
这告诉我的是整个矩阵看起来像这样
( 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 )
这意味着我可以只用零创建整个矩阵,然后更改对角线和相邻对角线。
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.
注意:我不熟悉这背后的数学原理,但J(n,m-1)
的值对我来说似乎是一个常数。
注2:我把值放在i + 1
索引处,因为我是从位置1;1
开始的,但你可以从相反的方向开始,i-1
。