2
$\begingroup$

I would like to construct a matrix $G$ composed of block matrices $G^{(i)}$ as defined entry-wise below stacked vertically as $$G = \begin{bmatrix} G^{(1)} \\ G^{(2)} \\ \vdots \\ G^{(i)} \\ \vdots \\ G^{(n)} \end{bmatrix}$$ for some natural number $n$.

$G^{(i)}_{j,k}$ stands for the element of the $j$'th row and the $k$'th column of the $i$'th block matrix. \begin{align} G^{(i)}_{j,k} :=& \frac{t_{i+1}-t}{h_i}\delta_{k,i}+\frac{t-t_i}{h_i}\delta_{k,i+1} -\frac16(t-t_i)(t_{i+1}-t)\Big(1+\frac{t_{i+1}-t}{h_i}\Big)\delta_{k,n+i} \\ &-\frac16(t-t_i)(t_{i+1}-t)\Big(1+\frac{t-t_i}{h_i}\Big)\delta_{k,n+i+1} \ \Bigg|_{t=t_{j,i}} \end{align} where $i\in\{1,2,\cdots,n\}$, $j\in\{1,2,\cdots,m(j)\}$, $m(j)$ is a natural numbered function of $j$, and $k\in\{1,2,\cdots,2n-2\}$, $\delta_{i,j}$ is the Kronecker delta function, and the vertical bar at the right end indicates $t$ is to be evaluated at $t_{j,t}$.

What is the most elegant and convenient way to do this?

$\endgroup$
  • $\begingroup$ ` t is to be evaluated at t_j,t` - explain please. You have a mixture of $t,t_i,t_{i+1}$ in your problem, but is $t$ a function? And if so why does it take either 1 (i or i+1) or two arguments (j,i)? $\endgroup$ – flinty Aug 4 at 17:49
  • $\begingroup$ @flinty: It just means substituting $t$ by $t_{j,i}$. $\endgroup$ – Hans Aug 4 at 19:23
  • $\begingroup$ Use SparseArray. Or DiagonalMatrix. What are the limits on j and k? Are the G's square matrices? Can n be anything or is it related to the bounds on i and j? More details would be helpful. $\endgroup$ – march Aug 4 at 21:08
  • $\begingroup$ You have a $\delta_{k,i}$ and a $\delta_{n+k,n+i}$, but $\delta_{n+k,n+i}=\delta_{k,i}$, so these actually are part of the same matrix element. I suspect that you meant something else, however, since you wrote them separately. Can you clarify? I think editing your post with a small example output would be very helpful, with explicit (small) limits on j and j like 3 or something, i = 2 (so, two blocks), and n = something. As it stands, $n$ actually doesn't do anything in this expression, as far as I can tell. $\endgroup$ – march Aug 4 at 21:19
  • $\begingroup$ @march: You are absolutely right. I have edited my question. Please review. Thank you. $\endgroup$ – Hans Aug 4 at 23:18
3
$\begingroup$

Assuming you have definitions already for n, single argument t[x], and two-argument t[x,y], then this is fairly straightforward:

g[i_, n_] := 
 Table[(t[i + 1] - t[j, i])/h[i] KroneckerDelta[k, 
     i] + (t[j, i] - t[i])/h[i] KroneckerDelta[k, i + 1] - 
   1/6 (t[j, i] - t[i]) (t[i + 1] - 
      t[j, i]) (1 + (t[i + 1] - t[j, i])/h[i]) KroneckerDelta[k, 
     n + i] - 
   1/6 (t[j, i] - t[i]) (t[i + 1] - 
      t[j, i]) (1 + ((t[j, i] - t[i])/h[i])) KroneckerDelta[k, 
     n + i + 1], {j, 1, 3}, {k, 1, 2*n-2}]

result = With[{n = 4}, Join @@ Array[g[#, n] &, n]]
result // Dimensions (*expected {12,3}*)
result // MatrixForm
| improve this answer | |
$\endgroup$
  • $\begingroup$ Is there a command to show the final explicit matrix form without the explicit KroneckerDelta function? $\endgroup$ – Hans Aug 4 at 19:32
  • 1
    $\begingroup$ @Hans can you explain what you're after? MatrixForm is just a presentational thing - you shouldn't use it in computations. $\endgroup$ – flinty Aug 4 at 22:05
  • $\begingroup$ Your solution is good. I understand I can use it for computation. I am wondering if there is a command to see the matrix in the explicit more conventional form as a sanity check. $\endgroup$ – Hans Aug 4 at 22:55
  • $\begingroup$ Yes that's what the //MatrixForm is for. It just looks weird because the expressions in the matrix are very wide. $\endgroup$ – flinty Aug 4 at 22:57
  • 1
    $\begingroup$ @Hans Array[g[#, n] &, n] generates a list of {g[1,n],g[2,n],g[3,n],..., g[n,n]} and @@ splices this list into a sequence of arguments for Join , so Join[g[1,n],g[2,n],g[3,n]...]. This is stacking them 'vertically' though Mathematica isn't like matlab; it works on lists, so it doesn't have a notion of vertical/horizontal, just lists and list nesting. If you wanted them horizontal you would use Join[g[1,n],g[2,n],g[3,n]....,{2}] with a level spec {2} at the end. $\endgroup$ – flinty Aug 5 at 11:26
2
$\begingroup$

Is this a formula with the Einstein sum convention. If so then the summation over two same indices has to be carried through.

g[i_, j_, k_, n_] := 
 Sum[(t[i + 1] - t[j, i])/h[ii] KroneckerDelta[k, 
     ii] + (t[j, i] - t[ii])/h[ii] KroneckerDelta[k, ii + 1] - 
   1/6 (t[j, i] - t[ii]) (t[ii + 1] - 
      t[j, i]) (1 + (t[ii + 1] - t[j, i])/h[ii]) KroneckerDelta[n + k,
      n + ii] - 
   1/6 (t[j, i] - t[ii]) (t[ii + 1] - 
      t[j, i]) (1 + ((t[j, i] - t[ii])/h[ii])) KroneckerDelta[n + k, 
     n + ii + 1], {ii, 1, 3}]

g[1, j, k, n]

(KroneckerDelta[1, k] (t[2] - t[j, 1]))/h[1] + (
 KroneckerDelta[2, k] (t[2] - t[j, 1]))/h[2] + (
 KroneckerDelta[3, k] (t[2] - t[j, 1]))/h[3] + (
 KroneckerDelta[2, k] (-t[1] + t[j, 1]))/h[1] - 
 1/6 KroneckerDelta[1 + n, 
   k + n] (1 + (t[2] - t[j, 1])/h[1]) (t[2] - t[j, 1]) (-t[1] + 
    t[j, 1]) + (KroneckerDelta[3, k] (-t[2] + t[j, 1]))/h[2] - 
 1/6 KroneckerDelta[2 + n, 
   k + n] (1 + (t[3] - t[j, 1])/h[2]) (t[3] - t[j, 1]) (-t[2] + 
    t[j, 1]) + (KroneckerDelta[4, k] (-t[3] + t[j, 1]))/h[3] - 
 1/6 KroneckerDelta[3 + n, 
   k + n] (1 + (t[4] - t[j, 1])/h[3]) (t[4] - t[j, 1]) (-t[3] + 
    t[j, 1]) - 
 1/6 KroneckerDelta[2 + n, 
   k + n] (t[2] - t[j, 1]) (-t[1] + t[j, 1]) (1 + (-t[1] + t[j, 1])/
    h[1]) - 1/
  6 KroneckerDelta[3 + n, 
   k + n] (t[3] - t[j, 1]) (-t[2] + t[j, 1]) (1 + (-t[2] + t[j, 1])/
    h[2]) - 1/
  6 KroneckerDelta[4 + n, 
   k + n] (t[4] - t[j, 1]) (-t[3] + t[j, 1]) (1 + (-t[3] + t[j, 1])/
    h[3])

The interpretation distinguishes between the indices i of the Einstein summation and the indices i from the discretization of time i in t_ij.

Example for displaying this as a matrix is

Table[h[i, j, k, n], {i, 1, 3}] // MatrixForm

enter image description here

This is much to lengthy for g.

| improve this answer | |
$\endgroup$
  • $\begingroup$ Actually, there is no summation over 'ii' here. 'ii' should be the index of the block matrix. $\endgroup$ – Hans Aug 4 at 19:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.