1
$\begingroup$

I am solving the matrix equation to pass from the Fokker-Planck to the Langevin equation. Specifically: \begin{equation} D=BB^T \end{equation} where $$ D= \left(\begin{matrix} l x_1(1-r_1) + e_1 r_1 & 0 & l x_1(1-r_1) & 0 \\ 0 & l x_2(1-r_2) + e_2 r_2 & 0 & l x_2(1-r_2) \\ l x_1(1-r_1) & 0 & l x_1+2D_{11}x_1 + D_{12}x_1 +D_{21}x_2 c_{11}r_1 +c_{12}r_2 & D_{12}x_1 +D_{21}x_2 \\ 0 & l x_2(1-r_2) & D_{12}x_1 +D_{21}x_2 & l x_2+2D_{22}x_2 + D_{12}x_1 +D_{21}x_2 c_{21}r_1 +c_{22}r_2 \end{matrix}\right) $$

where $B$ is symmetric. I am trying to solve the equation for $\mathbb{B}$ as follows:

Dmatrix = {{e1 r1 + l x1 (1 - r1), 0, l x1 (1 - r1), 0}, {0, e2 r2 + l x2 (1 - r2), 0, l x2 (1 - r2)}, {l x1 (1 - r1), 0, l x1 + 2 D11 x1 + c11 r1 + D12 x1 + D21 x2 + c12 r2, D12 x1 + D21 x2}, {0, l x2 (1 - r2), D12 x1 + D21 x2, l x2 + D21 x2 + 2 D22 x2 + D12 x1 + c21 r1 + c22 r2}}; B = Array[Y, {4, 4}]; For[i = 1, i < 5, i++, For[j = i, j < 5, j++, B[[j, i]] = B[[i, j]]]]; sol = Solve[Dmatrix == B . Transpose[B], Flatten[B]]; B /. Flatten[sol] // MatrixForm

But I am not finding a proper solution. Mathematically I know that more than a solution will be possible, since these are 16 quadratic equations. Can anyone help me?

Improve this question
$\endgroup$
2
  • $\begingroup$ This is a polynomial system of 16 equations in 16 unknowns with several parameters. I think no chance to find its symbolic solution. $\endgroup$ Commented Sep 19, 2023 at 16:12
  • 1
    $\begingroup$ Finding (what is effectively) the square root of a matrix, is a notoriously difficult problem. One may use a particular decomposition, eg Cholesky etc (see the help). $\endgroup$ Commented Sep 19, 2023 at 17:07

2 Answers 2

Reset to default
4
$\begingroup$

Doing a Cholesky decomposition, gives the answer practically instantaneously:

Dmatrix = {{e1 r1 + l x1 (1 - r1), 0, l x1 (1 - r1), 0}, {0, e2 r2 + l x2 (1 - r2), 0, l x2 (1 - r2)}, {l x1 (1 - r1), 0, l x1 + 2 D11 x1 + c11 r1 + D12 x1 + D21 x2 + c12 r2,D12 x1 + D21 x2}, {0, l x2 (1 - r2), D12 x1 + D21 x2,l x2 + D21 x2 + 2 D22 x2 + D12 x1 + c21 r1 + c22 r2}}; B = CholeskyDecomposition[Dmatrix] //. Conjugate[X_] -> X // Simplify (* {{Sqrt[e1 r1 + l x1 - l r1 x1], 0, -((l (-1 + r1) x1)/Sqrt[e1 r1 - l (-1 + r1) x1]), 0}, {0, Sqrt[ e2 r2 + l x2 - l r2 x2], 0, -((l (-1 + r2) x2)/Sqrt[e2 r2 - l (-1 + r2) x2])}, {0, 0, Sqrt[ c11 r1 + c12 r2 + 2 D11 x1 + D12 x1 + l x1 + ( l^2 (-1 + r1)^2 x1^2)/(-e1 r1 + l (-1 + r1) x1) + D21 x2], ( D12 x1 + D21 x2)/Sqrt[ c11 r1 + c12 r2 + 2 D11 x1 + D12 x1 + l x1 + ( l^2 (-1 + r1)^2 x1^2)/(-e1 r1 + l (-1 + r1) x1) + D21 x2]}, {0, 0, 0, \[Sqrt](c21 r1 + c22 r2 + D12 x1 + D21 x2 + 2 D22 x2 + l x2 - (D12 x1 + D21 x2)^2/( c11 r1 + c12 r2 + 2 D11 x1 + D12 x1 + l x1 + ( l^2 (-1 + r1)^2 x1^2)/(-e1 r1 + l (-1 + r1) x1) + D21 x2) + ( l^2 (-1 + r2)^2 x2^2)/(-e2 r2 + l (-1 + r2) x2))}} *)

The result looks a bit ugly, but it works:

Transpose[B].B - Dmatrix // Simplify (* {{0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}} *)

In the above I have made the very bold assumption that everything is real (note the Conjugate[X_] -> X part, so caveat emptor!

Finally, as mentioned previously this decomposition is not unique, however is extremely fast (0.00081 secs as reported by RepeatedTiming on my machine).

Improve this answer
$\endgroup$
0
1
$\begingroup$

This is only half the answer. I am able to solve for D, but to get B from D seems another problem.

Note, B is symmetric. But then B.Transpose is B.B. And this is also symmetric.

We may solve Dmatrix == D, where D is symmetric:

D = Array[Y, {4, 4}]; For[i = 1, i < 5, i++, For[j = i, j < 5, j++, D[[j, i]] = D[[i, j]]]]; sol = Solve[Dmatrix == B, DeleteDuplicates[Flatten[B]]][[1]]; B2 = B /. sol {{e1 r1 + l x1 - l r1 x1, 0, -l (-1 + r1) x1, 0}, {0, e2 r2 + l x2 - l r2 x2, 0, -l (-1 + r2) x2}, {-l (-1 + r1) x1, 0, c11 r1 + c12 r2 + 2 D11 x1 + D12 x1 + l x1 + D21 x2, D12 x1 + D21 x2}, {0, -l (-1 + r2) x2, D12 x1 + D21 x2, c21 r1 + c22 r2 + D12 x1 + D21 x2 + 2 D22 x2 + l x2}}

But now we need to solve B2 == B.B. For this MMA takes longer than I am willing to wait.

Improve this answer
$\endgroup$
4
  • 1
    $\begingroup$ After a while B=MatrixPower[B2, 1/2] evaluates... $\endgroup$ Commented Sep 19, 2023 at 17:02
  • 2
    $\begingroup$ @UlrichNeumann: Can you kindly present the result and timing? I ask it because that command crashes my comp in short time. $\endgroup$ Commented Sep 19, 2023 at 18:43
  • $\begingroup$ @user64494 Evaluation B=MatrixPower[B2, 1/2]; takes 25s on Mathematica v12.2, but I can't display the huge result in finite time $\endgroup$ Commented Sep 20, 2023 at 8:45
  • $\begingroup$ @user64494 B//LeafCount gives 182264749, probably to much to display without crash? $\endgroup$ Commented Sep 20, 2023 at 9:15

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.