Message Boards Message Boards

GROUPS:

[✓] Solve a system of coupled ODE's?

Posted 4 months ago
360 Views
|
2 Replies
|
2 Total Likes
|

Hello!

I am currently working on neutrino's oscillations and I need to solve a system of coupled ODE's. The problem is that I have never coded in Mathematica so my code is probably awful and... doesn't work at all. The system to solve is the following:enter image description here (at the end, one must take the modulus of the amplitude)

with everything a constant except V_e a function of time. This is my code which does not return the right solution (see below)

Gf = 0;

Xmax = 100000;

O13 = N[Degree, 8.54];
O23 = N[Degree, 47.2];
O12 = N[Degree, 33.62];
dcp = N[Degree, 234];

dm21 = 7.40*10^(-5);
dm31 = 2.494*10^(-3);

p = 1;

M = {{0, 0, 0}, {0, dm21, 0}, {0, 0, dm31}};

R1 = {{1, 0, 0}, {0, Cos[O23], Sin[O23]}, {0, -Sin[O23], Cos[O23]}};
R2 = {{Cos[O13], 0, Sin[O13]*Exp[-I*dcp]}, {0, 1, 
    0}, {-Sin[O13]*Exp[I*dcp], 0, Cos[O13]}};
R3 = {{Cos[O12], Sin[O12], 0}, {-Sin[O12], Cos[O12] , 0}, {0, 0, 1}};
U = R1*R2*R3;
Ve[x_] = Sqrt[2]*Gf*Exp[-2*x/Xmax];
P[x_] = 1/(2 p) * U*M*
    ConjugateTranspose[U] + {{Ve[x], 0, 0}, {0, 0, 0}, {0, 0, 0}};
A[x_] = Table[ai[i, x], {i, 1, 3}];
A0 = Table[0, {3}];
A0[[1]] = 1;

diffEq = Table[
  A'[x][[i]] == -I (P[x].A[x])[[i]], {i, 1, 3}]; initialCond = 
 Table[A[0][[i]] == A0[[i]], {i, 1, 3}];
sol = NDSolve[{Join[diffEq, initialCond]}, A[x], {x, 0, Xmax}];
Plot[Evaluate[
  Table[(ai[i, x] /. sol)*Conjugate[ai[i, x] /. sol], {i, 1, 3}]], {x,
   0, Xmax}]

The solution is supposed to be enter image description here

Thank you very much for your help !

2 Replies

I don't have time to track down the problem (at first glance, I'm not sure there's enough information given), but I can show you what ODE system you're integrating, and maybe you can track down the problem from there:

Join[diffEq, initialCond]
(*
{Derivative[0, 1][ai][1, x] == 0, 
 Derivative[0, 1][ai][2, x] == (0. - 0.00003697746403230496*I)* ai[2, x], 
 Derivative[0, 1][ai][3, x] == (0. - 0.0012462404769806562*I) * ai[3, x],
 ai[1, 0] == 1, ai[2, 0] == 0, ai[3, 0] == 0}
*)

From the first ODE, we can see that ai[1, x] should be a constant, and from the IC,that the constant is 1. The second and third ODEs are homogeneous linear, so the ICs imply the solutions ai[2, x] and ai[3, x] are the constant zero (or "trivial") solution. The system is a direct product of three independent ODEs. The desired solution looks like they should be coupled.

One thing that immediately jumps out are these assignments:

O13 = N[Degree, 8.54];
O23 = N[Degree, 47.2];
O12 = N[Degree, 33.62];
dcp = N[Degree, 234];

You might want to examine the actual values computed. My guess is that you wanted to multiply the numbers by Degree, instead of assigning all four parameters the same numeric value (with different precisions). For example, O13 = 8.54 * Degree.

Posted 4 months ago

Problem solved, I just didn't use the right operator for matrix multiplication :/

Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract