Message Boards Message Boards

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

GROUPS:

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 !

POSTED BY: Loic Sablon
Answer
1 month ago

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 BY: Michael Rogers
Answer
1 month ago

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

POSTED BY: Loic Sablon
Answer
1 month ago

Group Abstract Group Abstract