# Solving system of ODEs with extra parameter

I would like to solve a $$2\times 2$$ system of the form $$\frac{d}{d\theta}T=TA,\quad T(0)=Id$$ where $$\theta$$ is real and $$A$$ is of the form $$A=\begin{pmatrix} 0 & \frac{e^{-i \theta}}{\lambda}\\ \frac{1}{36}e^{-i\theta}\left(9\lambda + 2(\lambda-1)^2 (6\cos{\theta} + \cos{2\theta} + 6)\right) & 0\end{pmatrix},$$ with $$\lambda$$ a free parameter in the unit circle.

In particular I'm interested in obtaining numeric solutions at $$\theta=2\pi$$ depending on the extra parameter $$\lambda$$. I'm fairly new using Mathematica, and this is what I have tried so far:

T[θ_] = {{T11[θ], T12[θ]}, {T21[θ], T22[θ]}};
A[θ_] = {
{0, E^(-I θ)/λ},
{1/36 E^(-I θ) (9 λ + 2 (-1 + λ)^2 (6 + 6 Cos[θ] + Cos[2 θ])), 0}
};
sys = {T'[θ] == T[θ].A[θ]};


The previous code sets the system that I want to solve and now I try to solve numerically. First I've tried

NSol = NDSolve[
{sys, T11[0] == 1, T12[0] == 0, T21[0] == 0, T22[0] == 1},
{T11[θ], T12[θ], T21[θ], T22[θ]},
{θ},
{θ, 0, 2 Pi}
];


which gives me the output

NDSolve::dupv: "Duplicate variable θ found in NDSolve[<<1>>]."


I have also tried

Nsol2 = ParametricNDSolve[
{sys, T11[0] == 1, T12[0] == 0, T21[0] == 0, T22[0] == 1},
{T11, T12, T21, T22},
{θ, 0, 2 Pi},
{λ}
];


which gives me as output $$T_{11},\dots,T_{22}$$ as ParametricFunctions depending on each other and on $$\lambda$$.

I don't know if this is the right approach and, if so, how to extract a numeric expression depending on $$\lambda$$ from the last output - all that I've seen in the documentation are examples that are plotted for specific values of the parameter. Any help is much appreciated.

EDIT

Following the comments in one of the answers I'd like to explain further: the output that I would like to obtain is some sort of function depending of the parameter $$\lambda$$ that I can manipulate afterwards. Say for instance, computing the series expansion of powers of $$\lambda$$ of my solution. I don't know how to treat the ParametricFunction that I obtain to do such computations.

Here is how to obtain the series expansion of the matrix components. First, using a tweaked version (solving for t[2π] instead of t) of JM's formulation:

A[θ_] = {
{0, E^(-I θ)/λ},
{1/36 E^(-I θ) (9 λ + 2 (-1 + λ)^2 (6 + 6 Cos[θ] + Cos[2 θ])), 0}
};

pf = ParametricNDSolveValue[
{
t'[θ] == t[θ].A[θ], t[0] == IdentityMatrix[2]
},
t[2π],
{θ, 0, 2π},
λ
];


Then, pf will return the matrix value at $$\theta = 2 \pi$$. For example:

pf[1]
pf[Exp[I Pi/3]]


{{1. + 3.62818*10^-9 I, 6.82646*10^-8 - 5.73536*10^-9 I}, {1.70661*10^-8 - 1.43384*10^-9 I, 1. + 3.62818*10^-9 I}}

{{0.985595 + 1.17074 I, -0.425572 + 0.737112 I}, {-0.788363 - 1.36549 I, 0.985595 - 1.17074 I}}

Finding the series expansion is simple:

DecimalForm[Series[pf[λ], {λ, 1, 5}], {4,4}] //TeXForm


$$\left( \begin{array}{cc} 1.0000 & 6.2830 \\ 1.5710 & 1.0000 \\ \end{array} \right)+\left( \begin{array}{cc} 0.0000+0.0000 i & -6.2830+0.0000 i \\ 1.5710+0.0000 i & 0.0000+0.0000 i \\ \end{array} \right) (\lambda -1)+\left( \begin{array}{cc} 0.0000-0.7869 i & 0.6691-0.0000 i \\ 0.8799+0.0000 i & -0.0000+0.7869 i \\ \end{array} \right) (\lambda -1)^2+\left( \begin{array}{cc} 0.0000+0.7869 i & -1.3380+0.0000 i \\ 0.0000+0.0000 i & 0.0000-0.7869 i \\ \end{array} \right) (\lambda -1)^3+\left( \begin{array}{cc} -0.0152-0.4524 i & 1.8410-0.0000 i \\ -0.5708-0.0000 i & -0.0152+0.4524 i \\ \end{array} \right) (\lambda -1)^4+\left( \begin{array}{cc} 0.0305+0.1179 i & -2.1780+0.0000 i \\ 0.5708-0.0000 i & 0.0305-0.1179 i \\ \end{array} \right) (\lambda -1)^5+O\left((\lambda -1)^6\right)$$

As far as I can tell, you did everything right in your second approach. Classically, Mathematica returns lists of rules as results of solving functions, but in this case, I find this tradition rather confusing and prefer to use ParametricNDSolveValue; it returns a ParametricFunction object that, when applied to a numerical parameter, returns a list of 4 InterpolatingFunction for your 4 functions.

T[θ_] = {{T11[θ], T12[θ]}, {T21[θ], T22[θ]}};
A[θ_] = {{0, E^(-I θ)/λ}, {1/36 E^(-I θ) (9 λ + 2 (-1 + λ)^2 (6 + 6 Cos[θ] + Cos[2 θ])), 0}};
sys = {T'[θ] == T[θ].A[θ]};

Tsol = ParametricNDSolveValue[{sys, T11[0] == 1, T12[0] == 0,
T21[0] == 0, T22[0] == 1},
{T11, T12, T21, T22},
{θ, 0, 2 Pi},
{λ}
];


In order to obtain the numerical values for all the solutions at θ = 2 Pi for a given parameter, say λ = 0.1, you may use Through:

Through[Tsol[0.1][2. Pi]]


{-0.545795 + 1.00532 I, -1.43497 - 7.95125*10^-7 I, -0.215035 - 2.80298*10^-8 I, -0.545795 - 1.00532 I}

In order to make that into a function, you may use

f = λ \[Function] Through[Tsol[λ][2. Pi]]

• Thanks for your answer. Only the last step differs from what I'd like to obtain. What I would like to see is a function depending on $\lambda$, not the evaluation at a single point. – Edu Mar 30 '19 at 14:17

For this case, you don't even need to write out the components of your matrix function:

pf = ParametricNDSolveValue[{t'[θ] == t[θ].{{0, Exp[-I θ]/λ},
{Exp[-I θ] (9 λ + 2 (λ - 1)^2
(6 Cos[θ] + Cos[2 θ] + 6))/36, 0}},
t[0] == IdentityMatrix[2]}, t, {θ, 0, 2 π}, λ,
Method -> "StiffnessSwitching"];

sol = pf[(3 + 4 I)/5];

ParametricPlot[ReIm[Tr[sol[θ]]], {θ, 0, 2 π}]


ParametricPlot[ReIm[Det[sol[θ]]], {θ, 0, 2 π}]


You can even make a plot where the parameter is varying:

Plot[Re[Tr[pf[Exp[I ϕ]][2 π]]], {ϕ, 0, 2 π}]


• Thanks for your answer. Could I see the solution as a function of $\lambda$ instead? I suppose is something possible to do... – Edu Mar 30 '19 at 14:23
• That is what the third plot is; I let $\lambda=\exp(i\varphi)$ (quote "with $\lambda$ a free parameter in the unit circle") and plotted the real part of the trace of the matrix. – J. M.'s technical difficulties Mar 30 '19 at 14:24
• Sure, but I'd like to manipulate it further. Not just see its plot. See my point? – Edu Mar 30 '19 at 14:26
• You make no mention of what kind of manipulations you like to do in your question, so no, I do not see. – J. M.'s technical difficulties Mar 30 '19 at 14:26
• For example, can I get the solution as a function of $\lambda$? Or as a series expansion of powers of $\lambda$? – Edu Mar 30 '19 at 14:41