second order direction field in matlab

Question

I am trying to plot a second order direction field in MATLAB. The equation is $$LQ'' + RQ' + (1/C)Q = E(t)$$ where Q is the total charge on the capacitor at time t in coulombs, L is inductance in henrys, R is resistance in ohms, and E is impressed voltage in volts. If possible, I want to plot the particular solution, $$Q'' + 40Q' + 625Q = 100\cos(10t)$$ on top of the direction field. I don't really understand MATLAB, so I really appreciate any help anyone can offer here.

Answer 1

In order to get a plottable solution your first need to specify initial conditions for $Q$ and $\dot{Q}$. I assumed them to be $Q_0=1$ and $\dot{Q}_0=0$.

First rewrite the ODE as a first order system. Introduce the substitution $Q = q_1$ and $q_2 = \dot{q}_1$. This will transform the equation into:

$$ \dot{q}_1 = q_2$$ $$ \dot{q}_2 = -625q_1-40q_2+100\cos(10t).$$

Now create a function in MATLAB and save it as electric.m in your current folder:

function Dq = electric(t,q)    
%q(1),q(2) are the states
%Dq(i) is the derivative of the ith state
Dq = zeros(2,1);
Dq(1) = q(2);
Dq(2) = -625*q(1)-40*q(2)+100*cos(10*t);
end 

You should now be able to run the following lines from your command window:

% [0,10] time interval; [1;0] initial conditions for q1 and q2
[t,q] = ode45(@electric,[0,10],[1;0]); 
plot(t,q(:,1)); % plot q1=Q