I would like to know how to solve this system of differential equation.
It consist of 3 ODEs, describing the behavior of an Induction Machine supplied with DC voltage. I a interested to derive the input current to the machine(I1) but not numerically. In Matlab i could find the numerical solution and could plot the curve(Transferring to Laplace, solving equation and then doing inverse Laplace), but i need the answer to be a mathematics expression. I believe the answer should contain an exponential term+some hyperbolic term.
Laplace transform:(Maybe i don't really need to solve it through conversion to Laplace, because finding the inverse Laplace is very difficult)
Matlab output waveform:
Circuit:
clear;
syms r1 l1 r2 l2 rc lm real; syms i1(t) i2(t) i3(t) s; di1(t) = diff(i1(t),t); di2(t) = diff(i2(t),t); di3(t) = diff(i3(t),t);
E = 1; % Voltage eq1(t) = di1(t) + ((rc+r1)/l1)*i1(t) - (rc/l1)*i2(t)-E/l1; eq2(t) = di2(t) - di3(t) + (rc/lm)*i2(t) - (rc/lm)i1(t); eq3(t) = di3(t) - (lm/(lm+l2)) di2(t) + (r2/(lm+l2))*i3(t);
L1(t) = laplace(eq1,t,s); L2(t) = laplace(eq2,t,s); L3(t) = laplace(eq3,t,s); %{0.3,0.2,100,0.1,0.00386,0.004,0,0,0});
syms Li1 Li2 Li3 % Ni1=L1(t); % Ni2=L2(t); % Ni3=L3(t);
Ni1 = subs(L1(t),{r1,r2,rc,lm,l1,l2,i1(0),i2(0),i3(0)}, ... {0.3,0.2,100,0.1,0.00386,0.004,0,0,0});
Ni2 = subs(L2(t),{r1,r2,rc,lm,l1,l2,i1(0),i2(0),i3(0)}, ... {0.3,0.2,100,0.1,0.00386,0.004,0,0,0});
Ni3 = subs(L3(t),{r1,r2,rc,lm,l1,l2,i1(0),i2(0),i3(0)}, ... {0.3,0.2,100,0.1,0.00386,0.004,0,0,0});
Ni1 =... subs(Ni1,{laplace(i1(t),t,s),laplace(i2(t),t,s),laplace(i3(t),t,s)},{Li1,Li2,Li3}); NI1 = collect(Ni1,Li1)
Ni2 = ... subs(Ni2,{laplace(i1(t),t,s), laplace(i2(t),t,s),laplace(i3(t),t,s)}, {Li1,Li2,Li3});
Ni2 = collect(Ni2,Li2);
Ni3 = ... subs(Ni3,{laplace(i1(t),t,s), laplace(i2(t),t,s),laplace(i3(t),t,s)}, {Li1,Li2,Li3});
Ni3 = collect(Ni3,Li2);
[Li1, Li2, Li3] = solve(Ni1, Ni2, Ni3, Li1, Li2, Li3)
i1 = ilaplace(Li1, s, t); i2 = ilaplace(Li2, s, t); i3 = ilaplace(Li3, s, t);
simplify(i1); simplify(i2); simplify(i3);
i1=vpa(i1) i2=vpa(i2) i3=vpa(i3)
Di1=diff(i1); Ll1=E./Di1; Lxpll2=(diff(Di1));
t=0:1e-4:10; Ii1=subs(i1,t); DDi1=subs(Di1,t); LLl1=subs(Ll1,t); LLxpll2=subs(Lxpll2,t);
% i1=vpa(i1) % i2=vpa(i2)
% % % subplot(2,2,1); plot(t,Ii1); title('Input current'); ylabel('i1(t)[A]'); xlabel('t[s]');grid subplot(2,2,2); plot(t,DDi1); title('Derivation of the current: di1(t)/d(t)'); ylabel('di1(t)/d(t)');xlabel('t[s]'); grid subplot(2,2,3); plot(t,LLl1); title('Input inductance'); ylabel('Li[H]');xlabel('t[s]'); grid subplot(2,2,4); plot(t,LLxpll2); title('Second derivation of the current: d^2i1(t)/d(t)'); ylabel('d^2i1(t)/d(t)');xlabel('t[s]'); grid
Thanks in advance,
You have a series circuit of three things to which a signal is applied. $$ R1+L1+(RC || Lm || (R2+L2)), $$ where I am using $||$ to indicate parallel. This one circuit gives the main equation. Then you figure out how to look at any one of the three series items in this through a divider, and then burrow down into the parallel currents if needed.
More detail: I'm using your Laplace transform equations: $$ \begin{array}{llll} (sL_1+R_1)I_1 & +R_c(I_1-I_2) & & = E/s \\ & +R_c(I_2-I_1) & +sL_m(I_2-I_3) & = 0 \\ & & +sL_m(I_3-I_2) +(sL_2+R_2)I_3 & = 0 \end{array} $$ Add these three equations to relate $I_1$ to $I_3$: $$ (sL_1+R_1)I_1 +(sL_2+R_2)I_3 = E/s,\\ I_3 = \frac{E/s-(sL_1+R_1)I_1}{sL_2+R_2}. $$ Use this to relate $I_2$ to $I_1$ using the second loop equation: $$ \begin{align} (R_c +sL_m) I_2 & = R_c I_1+sL_m I_3 \\ & =R_c I_1 + sL_m\frac{E/s-(sL_1+R_1)I_1}{sL_2+R_2}\\ & =\left[R_c-\frac{sL_m(sL_1+R_1)}{sL_2+R_2}\right]I_1+\frac{L_m}{sL_2+R_2}E \end{align} $$ Therefore, $$ I_2 = \frac{1}{sL_m+R_c}\left[\left[R_c-\frac{sL_m(sL_1+R_1)}{sL_2+R_2}\right]I_1+\frac{L_m}{sL_2+R_2}E\right]\\ = \frac{R_c(sL_2+R_2)-sL_m(sL_1+R_1)}{(sL_m+R_c)(sL_2+R_2)}I_1 +\frac{L_m}{(sL_m+R_c)(sL_2+R_2)}E\\ = -\frac{L_m L_1 s^{2}+(L_m R_1-L_2 R_c)s-R_c R_2}{(sL_m+R_c)(sL_2+R_2)}I_1 +\frac{L_m}{(sL_m+R_c)(sL_2+R_2)}E $$ Rewrite the first loop equation $$ (sL_1+R_1+R_c)I_1 = R_c I_2 +E/s \\ I_1 = \frac{R_c}{sL_1+R_1+R_c}I_2+\frac{E}{s(sL_1+R_1+R_c)} \\ = -\frac{L_m L_1 s^{2}+(L_m R_1-L_2 R_c)s-R_c R_2}{(sL_m+R_c)(sL_2+R_2)(sL_1+R_1+R_2)}R_cI_1\\ - \frac{1}{(sL_m+R_c)(sL_2+R_2)(sL_1+R_1+R_2)}R_cL_mE \\ + \frac{E}{s(sL_1+R_1+R_c)} $$ Now you have $I_1$ in terms of $E$. $$ \left[1+\frac{L_m L_1 s^{2}+(L_m R_1-L_2 R_c)s-R_c R_2}{(sL_m+R_c)(sL_2+R_2)(sL_1+R_1+R_2)}R_c\right]I_1 \\ = \left[-\frac{R_cL_m}{(sL_m+R_c)(sL_2+R_2)(sL_1+R_1+R_2)} +\frac{1}{s(sL_1+R_1+R_c)}\right]E $$ Multiplying both sides by $(sL_m+R_c)(sL_2+R_2)(sL_1+R_1+R_2)$, simplifying the cubic coefficient $C(s)$ of $I_1$ and dividing by that coefficient gives $I_1 = \frac{1}{C(s)}\frac{p(s)}{s}E$, where $p$ is a cubic function polynomial. There are explicit solutions for the roots of a cubic. So, you can explicitly find the inverse Laplace transform after factoring $C(s)$ by using the inverse Laplace transform integral and the calculus of residues.