The title is a little ambiguous, but I didn't know how else to put it. What I'm trying to do is solve a system of equations

system1 = 
  {-I*ω*a1 == 
     -I*ω1*a1 - I*J12*a2 - I*J13*a3 - κ1[ω]/2*a1 - γ1[ω]/2*a1 + Sqrt[κ1[ω]]*ain, 
   -I*ω*a2 == -I*ω2*a2 - I*J12*a1 - I*J23*a3 - γ2/2*a2, 
   -I*ω*a3 == -I*ω3*a3 - I*J13*a1 - I*J23*a2 - γ3/2*a3}; 

Which I then solve for a1, a2 and a3:

s = Solve[system1, {a1, a2, a3}];

And then, finally, I am interested in a function $R(\omega)$ which is given by

I*Sqrt[κ1]*(a1 /. s)/ain - 1

If I simply use all the above and evaluate that last term, it works. I get a (rather ugly) expression of the type

$\left\{-1+\frac{2 i \text{$\kappa $1} \left(4 \text{J23}^2+(\text{$\gamma $2}-2 i (\omega -\text{$\omega $2})) (\text{$\gamma $3}-2 i (\omega -\text{$\omega $3}))\right)}{4 \text{J12}^2 (\text{$\gamma $3}-2 i (\omega -\text{$\omega $3}))-16 i \text{J12} \text{J13} \text{J23}+4 \text{J13}^2 (\text{$\gamma $2}-2 i (\omega -\text{$\omega $2}))+(\text{$\gamma $1}+\text{$\kappa $1}-2 i \omega +2 i \text{$\omega $1}) \left(4 \text{J23}^2+(\text{$\gamma $2}-2 i \omega +2 i \text{$\omega $2}) (\text{$\gamma $3}-2 i \omega +2 i \text{$\omega $3})\right)}\right\}$

(I apologize for the format, I don't know how to give variables and parameters subscripts in mathematica so for example $J_{12}$ comes out as $J12$)

However, I want to define this as a function, and I am not sure if

R[ω_] := I*Sqrt[κ1]*(a1 /. s)/ain - 1

does the trick.

Now, the real problem comes from the fact that I don't want just this analytic expression, I want to plot $R(\omega)$. The real and imaginary parts, as well as the argument. For this I need values for my constants of course, which is fine. But the problems start appearing when I consider for example $\kappa_1$, which is not a constant but a function. To be clear, every term in these equations, apart from $a_1$, $a_2$, $a_3$,$\kappa_1$ and $\gamma_1$, are constants. These other five depend on the only variable in the equations, $\omega$.

The main issue is that I don't know how I should do this. I can define for example $\kappa_1$ the way I want to without any problems using

κ1[ω_] := 1/(ω*Zc*(C1 + CJ1 + CJ3));

This gives me a function $\kappa_1$ that I can evaluate. But how do I put this into my system of equations, system1? Do I have to replace κ1 with κ1[ω]?

If I do that, then I run into trouble when solving the system and defining an equation. If anyone coudld help, I'd be very grateful. For completeness, below is the full list of parameters:

C1 = 190*10^-15;
C2 = 240*10^-15;
C3  = 270*10^-15;
L1 = 1.74*10^ -9;
L2 = L1;
L3 = L1;
CJ1 = 30*10^-15;
CJ2 = CJ1;
CJ3 = CJ1;
R1 = 2*10^5;
R2 = R1;
R3 = R1;
Zc = 50;
Cmatrix = {{C1 + Cκ, CJ1, CJ3}, {CJ1, C2, CJ2}, {CJ3, CJ2, C3}};
Lmatrix = {{L1, 0, 0}, {0, L1, 0}, {0, 0, L1}};
Linv = Inverse[Lmatrix];
Cinv = Inverse[Cmatrix];
ω1 = Sqrt[Cinv[[1, 1]]*Linv[[1, 1]]];
ω2 = Sqrt[Cinv[[2, 2]]*Linv[[2, 2]]];
ω3 = Sqrt[Cinv[[3, 3]]*Linv[[3, 3]]];
J12 = 0.5*Cinv[[1, 2]]/(Sqrt[Cinv[[1, 1]]]*Cinv[[2, 2]])*Sqrt[ω1*ω2];
J13 = 0.5*Cinv[[1, 3]]/(Sqrt[Cinv[[1, 1]]]*Cinv[[3, 3]])*Sqrt[ω1*ω3];
J23 = 0.5*Cinv[[2, 3]]/(Sqrt[Cinv[[2, 2]]]*Cinv[[3, 3]])*Sqrt[ω3*ω2];
Cκ = 80*10^-15;
γ2 = 1/(R2*(C2 + CJ1 + CJ2));
γ3 = 1/(R3*(C3 + CJ2 + CJ3));
κtilde = 1/(Zc*Cκ);
BigK[ω_] := 1 + (κtilde/ω)^2;
Cκtilde[ω_] := Cκ*(1 - 1/BigK[ω]);
κ1[ω_] := 1/(BigK[ω]*Zc*(C1 + CJ1 + Cκtilde[ω] + CJ3));
γ1[ω_] := 1/(R1*(C1 + CJ1 + CJ3 + Cκtilde[ω]));
  • $\begingroup$ I'm not sure I fully understand your question, but perhaps R[ω_] = Abs[(a1 /. s)/ain - 1] is what you want. Note I used Set, not SetDelayed. $\endgroup$
    – m_goldberg
    Jun 2, 2015 at 12:06
  • $\begingroup$ You are right, it's too vague. I've updated my main post! $\endgroup$
    – user129412
    Jun 2, 2015 at 12:11
  • $\begingroup$ I still can't figure out what ain is. $\endgroup$
    – m_goldberg
    Jun 2, 2015 at 12:22
  • $\begingroup$ Right, a product of terrible notation. Its $a_{in}$, a constant value that gets divided out at the end so it doesn't need an actual value. It is in there all the way until I get to $R(\omega)$, then it disappears. Context wise, it has to do with an incoming wave in an electronic circuit containing inductors, resistors and capacitors. I'm calculating the reflection coefficient. $\endgroup$
    – user129412
    Jun 2, 2015 at 12:23

1 Answer 1


Using your definitions of system1 and s, when I write

R[ω_] = Simplify[I*Sqrt[κ1[ω]]*(a1 /. s)/ain - 1][[1]]

I get

R[ω_] = 
  -1 + (J13*(4*J23^2 + (γ2 - (2*I)*(ω - ω2))*(γ3 - (2*I)*(ω - ω3)))*κ1[ω])/
    (4*((J12*J23 + J13*((I/2)*γ2 + ω - ω2))*
      (-J13^2 + ((I*γ3 + 2*ω - 2*ω3)*(2*ω - 2*ω1 + I*γ1[ω] + I*κ1[ω]))/4) - 
        ((2*J13*J23 + J12*(I*γ3 + 2*ω - 2*ω3))*(2*J12*J13 + 2*J23*ω - 
          2*J23*ω1 + I*J23*γ1[ω] + I*J23*κ1[ω]))/4))

With this definition of R, Re[R[ω]] ,Im[R[ω]], and Arg[R[ω]] all seem eminently plot-able as functions of ω. If this doesn't work for you, you need to tell us what goes wrong.

Please note the use of Set ( = ). You don't want to use SetDelayed for defining R. In this case, the righthand side of the definition needs to be evaluated when the definition is evaluated.

  • $\begingroup$ Thanks for the answer. Indeed, with your comment on not using :=, I get to that point. I will actually do a test with a system for which I do know the solution to check if it is not some numerical anomaly instead. The := seems to have ruined a lot.. $\endgroup$
    – user129412
    Jun 2, 2015 at 13:45
  • $\begingroup$ Something is definitely still wrong, but I can't figure out what or where. In any case, I've used a different approach now, and that does work. Thanks for your time though, you've answered all that was possible in the current question! $\endgroup$
    – user129412
    Jun 2, 2015 at 20:05

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.