# Solving Coupled system of differential equation with NDSolve

I am solving a system with coupled system of differential equation to obtain the time evolution of two vectors. The problem stated here is equivalent to solving Landau Lifschitz Gilbert equationLLG Equation for two systems with a coupling between them. In my problem, m1[t] and m2[t] are two 3 component column vectors in cartesian coordinates. m1[t]={m1x,m1y,m1z} and m2[t]={m2x,m2y,m2z} respectively. m1[t] and m2[t] are acted upon by a term called effective field (HAnis).This HAnis is derived from an energy term

EAnis[r, \[Theta], \[CurlyPhi]] := (K1*(Sin[\[Theta]]^2) +
K2*(Sin[\[Theta]]^4) + K3*(Sin[\[Theta]]^6)*(Cos[6 \[CurlyPhi]]));
HAnis = Simplify[
r*Grad[EAnis[r, \[Theta], \[CurlyPhi]], {r, \[Theta], \[CurlyPhi]},
"Spherical"]]/MS;


Theta and phi are the polar and azimuthal angles in a spherical polar coordinate geometry. The time derivative of m1[t] and m2[t] is given by

cons1[t_] := -gamma*(10^(-9))*(-Hexch*Cross[m1[t], m2[t]] +
Cross[m1[t], (HAnis + Hext)]);
tGilbdamp1[t_] := alphag*Cross[m1[t], cons1[t]];
LLGS1 = {m1'[t] == cons1[t] + tGilbdamp1[t],
m1[0] == ToSphericalCoordinates[{1, 0, 0}]};


and

cons2[t_] := -gamma*(10^(-9))*(-Hexch*Cross[m2[t], m1[t]] +
Cross[m2[t], (HAnis + Hext)]);
tGilbdamp2[t_] := alphag*Cross[m2[t], cons2[t]];
LLGS2 = {m2'[t] == cons2[t] + tGilbdamp2[t],
m2[0] == ToSphericalCoordinates[{-1, 0, 0}]};


I tried to obtain numerical solution by

    sol1 = NDSolve[{LLGS1, LLGS2}, {m1, m2}, {t, tStart, tFinish},
StartingStepSize -> tStep,
Method -> {"FixedStep",
Method -> "ExplicitRungeKutta"},
AccuracyGoal -> 10, PrecisionGoal -> 100, MaxSteps -> \[Infinity]];


And I have stumbled upon the following error messages.

Cross::nonn1: The arguments are expected to be vectors of equal length, and the number of arguments is expected to be 1 less than their length.

Cross::nonn1: The arguments are expected to be vectors of equal length, and the number of arguments is expected to be 1 less than their length.

Cross::nonn1: The arguments are expected to be vectors of equal length, and the number of arguments is expected to be 1 less than their length.

General::stop: Further output of Cross::nonn1 will be suppressed during this calculation.

NDSolve::ndnum: Encountered non-numerical value for a derivative at t == 0.


Any help will be highly appreciated. The constants in the above equations are as follows:

gamma = 1.760859*(10^11);
alphag = (10^(-4));
tStep = (10^(-3));
tStart = 0;
tFinish = 8;
K1 = (3*(10^5));
K2 = (5*(10^5));
K3 = (5*(10^5));
Hext = {0, 0, 0};
Hexch = 1000;
MS = 238732.41;

• Please give all necessary information to run your code. For now, Hexch, Hext, and tSlondamp2[t] are undefined. Also, sometimes you use m1 and m2 without explicit time-dependence (m1[t] and m2[t]). – Roman Mar 12 '19 at 9:13
• Sorry @Roman. Thanks for your reply. tSlondamp[t] can be discarded for now and I have updated the code in the post. The values of Hexch (scalar) and Hext (3 component coloumn vector) are also updated. Hexch = 1000 and Hext = {0,0,0}. m1 and m2 are time dependent vectors. – Matcax Mar 12 '19 at 9:23
• I think you're mixing up Cartesian and spherical coordinates. ToSphericalCoordinates[{1, 0, 0}] gives spherical coordinates, but you equate it to m1[0] which you say is in Cartesian coordinates. Same for Hanis being in spherical coordinates but then appears in a Cross product with m1[t] (in the definition of cons1). As the code stands I think it cannot work, even if you fix the error messages. – Roman Mar 12 '19 at 9:43

If we correct all typos and add data for $$\theta, \phi$$, then the code works

gamma = 1.760859*10^11;
alphag = 10^(-4);
tStep = 10^(-3);
tStart = 0;
tFinish = 8;
K1 = 3*10^5;
K2 = 5*10^5;
K3 = 5*10^5;
Hext = {0, 0, 0};
Hexch = 1000;
MS = 238732.41; EAnis[r_, \[Theta]_, \[CurlyPhi]_] := K1*Sin[\[Theta]]^2 + K2*Sin[\[Theta]]^4 + K3*Sin[\[Theta]]^6*Cos[6*\[CurlyPhi]];
HAnis = Simplify[r*Grad[EAnis[r, \[Theta], \[CurlyPhi]], {r, \[Theta], \[CurlyPhi]}, "Spherical"]]/MS;
cons1[t_] := ((-gamma)*((-Hexch)*Cross[m1[t], m2[t]] + Cross[m1[t], HAnis + Hext]))/10^9;
tGilbdamp1[t_] := alphag*Cross[m1[t], cons1[t]];
LLGS1 = {Derivative[1][m1][t] == cons1[t] + tGilbdamp1[t], m1[0] == ToSphericalCoordinates[{1, 0, 0}]};
cons2[t_] := ((-gamma)*((-Hexch)*Cross[m2[t], m1[t]] + Cross[m2[t], HAnis + Hext]))/10^9;
tGilbdamp2[t_] := alphag*Cross[m2[t], cons2[t]];
LLGS2 = {Derivative[1][m2][t] == cons2[t] + tGilbdamp2[t], m2[0] == ToSphericalCoordinates[{-1, 0, 0}]};
sol1 = NDSolve[{LLGS1,
LLGS2} /. {\[Theta] -> 1, \[CurlyPhi] -> 0}, {m1, m2}, {t, tStart,
tFinish}]

Plot[m1[t] /. sol1, {t, tStart, tFinish}]

Plot[m2[t] /. sol1, {t, tStart, tFinish}]


• As I wrote above, this code mixed Cartesian and spherical coordinates erratically. Now it runs, yes; but the output is probably useless. – Roman Mar 12 '19 at 12:48
• @Roman We are fixing code errors here. Errors model let the author fixes himself. – Alex Trounev Mar 12 '19 at 12:52
• Yes Alex I agree. – Roman Mar 12 '19 at 13:54
• @Alex and Roman Thanks for your help. I have tried in a different way. I wrote m1 and m2 in cartesian co-ordinates. Used spherical polar derivatives for E and performed orthogonal basis transformation to convert them into cartesian coordinates. And then applied used ND Solve to obtain a solution. – Matcax Mar 15 '19 at 7:42
• @Matcax What does this have to do with the Landau-Lifschitz-Gilbert model? – Alex Trounev Mar 15 '19 at 14:14