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?

$\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$
    – user64494
    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$
    – Hans Olo
    Sep 19, 2023 at 17:07

2 Answers 2

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).

$\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.

$\endgroup$
4
  • 1
    $\begingroup$ After a while B=MatrixPower[B2, 1/2] evaluates... $\endgroup$ 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$
    – user64494
    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$ Sep 20, 2023 at 8:45
  • $\begingroup$ @user64494 B//LeafCount gives 182264749, probably to much to display without crash? $\endgroup$ Sep 20, 2023 at 9:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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