[ Contructing a triadiagonal matrix in Mathematica where nonzero elements contain functions and variables ]
Suppose I want to construct a matrix A
such that A[[i,i]]=f[x_,y_]+d[i]
, A[[i,i+1]]=u[i]
, A[[i+1,i]]=l[i], i=1,N
. Say, f[x_,y_]=x^2+y^2
.
How can I code this in Mathematica?
Additionally, if I want to integrate the first diagonal element of A
, i.e. A[[1,1]]
over x
and y
, both running from 0 to 1, how can I do that?
Answer 1
In[1]:= n = 4;
f[x_, y_] := x^2 + y^2;
A = Normal[SparseArray[{
{i_,i_}/;i>1 -> f[x,y]+ d[i],
{i_,j_}/;j-i==1 -> u[i],
{i_,j_}/;i-j==1 -> l[i-1],
{1, 1} -> Integrate[f[x,y]+d[1], {x,0,1}, {y,0,1}]},
{n, n}]]
Out[3]= {{2/3+d[1], l[1], 0, 0},
{u[1], x^2+y^2+ d[2], l[2], 0},
{0, u[2], x^2+y^2+d[3], l[3]},
{0, 0, u[3], x^2+y^2+d[4]}}
Answer 2
Band
is tailored specifically for this:
myTridiagonalMatrix@n_Integer?Positive :=
SparseArray[
{ Band@{1, 1} -> f[x, y] + Array[d, n]
, Band@{1, 2} -> Array[u, n - 1]
, Band@{2, 1} -> Array[l, n - 1]}
, {n, n}]
Check it out (no need to define f
, d
, u
, l
):
myTridiagonalMatrix@5 // MatrixForm
Note that MatrixForm
should not be part of a definition. For example, it's a bad idea to set A = (something) // MatrixForm
. You will get a MatrixForm
object instead of a table (= array of arrays) or a sparse array, and its only purpose is to be pretty-printed in FrontEnd. Trying to use MatrixForm
in calculations will yield errors and will lead to unnecessary confusion.
Integrating the element at {1, 1}
:
myTridiagonalMatrixWithFirstDiagonalElementIntegrated@n_Integer?Positive :=
MapAt[
Integrate[#, {x, 0, 1}, {y, 0, 1}]&
, myTridiagonalMatrix@n
, {1, 1}]
You may check it out without defining f
or d
, as well:
myTridiagonalMatrixWithFirstDiagonalElementIntegrated@5
The latter operation, however, looks suspicious. For example, it does not leave your matrix (or its corresponding linear system) invariant w.r.t. reasonable transformations. (This operation does not even preserve linearity of matrices.) You probably don't want to do it.
Comment on comment above: there's no need to define A[x_, y_] := …
to Integrate[A[[1,1]], {x,0,1}, {y,0,1}]
. Note that A[[1,1]]
is totally different from A[1, 1]
: the former is Part[A, 1, 1]
which is a certain element of table A
. A[1, 1]
is a different expression: if A
is some table then A[1, 1]
is (that table)[1, 1]
, which is a valid expression but is normally considered meaningless.