1
$\begingroup$

I am trying to input a difference equation in mathematica. The equation is $u_{x+1,y}+u_{x-1,y}+u_{x,y+1}+u_{x,y-1}+(\omega^2-4)u_{x,y}=\delta_{x,x0}\delta_{y,y0}$, where $u_{x,y}\in\mathbb{C}$,$(x,y)\in\mathbb{Z}^2$ is unknown and $\omega\in\mathbb{C}$, $x_0\in\mathbb{Z}$ and $y_0\in\mathbb{Z}$ are constants. This gives rise to a system of algebraic equations, $Ax=b$ type that can be solved. However, I don't have an efficient way to enter these equations in mathematica. I have used KroneckerDelta and SparseArray functions to make the matrices A and b; however, it takes a lot of time to form these matrices if the mesh size is large, say $(2n+1)\times(2n+1)$, where $n$ can be 200. Please help me out here. Here is my code

Amatrix[j_] := 
Table[KroneckerDelta[n, m] N[\[Omega]^2 - 4] + 
KroneckerDelta[n - 1, m] + KroneckerDelta[n, m - 1] + 
KroneckerDelta[n + j - 2, m] + KroneckerDelta[n, m + j - 2] + 
KroneckerDelta[n + (j - 2)^2 + 2*(j - 2) - (j - 2), m] + 
KroneckerDelta[n, m + (j - 2)^2 + 2*(j - 2) - (j - 2)], {n, 
1, (j - 2)^2 + 2*(j - 2)}, {m, 1, (j - 2)^2 + 2*(j - 2)}];

Bmatrix[j_] := 
Flatten@Table[
KroneckerDelta[x, x0] KroneckerDelta[y, y0], {y, 1, j}, {x, 1, 
 j - 2}];

Here, $j$ represents the mesh size. Here $2n+1$. Thank you.

$\endgroup$
  • $\begingroup$ Have you tried RSolve? $\endgroup$ – David Baghdasaryan May 15 '17 at 11:34
  • $\begingroup$ Yes, I did. It doesn't work in my case. $\endgroup$ – Gaurav Maurya May 15 '17 at 11:49
  • $\begingroup$ Can you please show the code that creates matrix of coefficients A and vector b $\endgroup$ – David Baghdasaryan May 15 '17 at 11:58
  • 1
    $\begingroup$ Hi, welcome to Mathematica.SE, please take the tour so you learn the basic rules of the site. Here its considered helpful to share your code attempts AND representative data in a well formatted form, so we can quickly see the problem you are facing. Examples of how to share your data here.Please edit your question and include the code of your own attempts in the question, and not in the comments $\endgroup$ – rhermans May 15 '17 at 12:31
  • $\begingroup$ You appear to wish to solve a finite difference approximation to the Helmholtz equation to obtain a numerical Green's function. If so, you also need boundary conditions. Perhaps, though, it would be better to use built-in functions, like GreenFunction, and then discretize the result. $\endgroup$ – bbgodfrey May 15 '17 at 13:51
1
$\begingroup$

simply forming the arrays, SparseArray is considerably faster.

n = 20
\[Omega] = .1
Amatrix[j_] := Table[
   KroneckerDelta[n, m] N[\[Omega]^2 - 4] +
    KroneckerDelta[n - 1, m] +
    KroneckerDelta[n, m - 1] +
    KroneckerDelta[n + j - 2, m] +
    KroneckerDelta[n, m + j - 2] +
    KroneckerDelta[n + (j - 2)^2 + 2*(j - 2) - (j - 2), m] + 
    KroneckerDelta[n, m + (j - 2)^2 + 2*(j - 2) - (j - 2)], {n, 
    1, (j - 2)^2 + 2*(j - 2)}, {m, 1, (j - 2)^2 + 2*(j - 2)}];

AmatrixS[j_] :=
 SparseArray[{
   {n_, m_} /; n == m -> N[\[Omega]^2 - 4],
   {n_, m_} /; n == m + 1 -> 1,
   {n_, m_} /; n == m - 1 -> 1,
   {n_, m_} /; n + j - 2 == m -> 1,
   {n_, m_} /; n == m + j - 2 -> 1,
   {n_, m_} /; n + (j - 2)^2 + 2*(j - 2) - (j - 2) == m -> 1,
   {n_, m_} /; n == m + (j - 2)^2 + 2*(j - 2) - (j - 2) -> 1
   }, {(j - 2)^2 + 2*(j - 2), (j - 2)^2 + 2*(j - 2)}]

sparse = AmatrixS[2 n + 1]; // AbsoluteTiming //First

0.0481058

original = Amatrix[2 n + 1] ; // AbsoluteTiming //First

35.6627

sparse == original

True

there is also considerable size savings as the array is extremely sparse:

ByteCount /@ {original, sparse}

{62220816, 272856}

notice I assumed your KroneckerDelta's are never nonzero at the same position. If that's not the case you would create a separate SparseArray for each delta and add them together.

$\endgroup$
  • $\begingroup$ Hey @george2079! This looks great. Thank you. I previously used SpareArray but not this way. Thank you again. $\endgroup$ – Gaurav Maurya May 17 '17 at 11:05

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.