User:Egm6321.f12.team4.mishra/Homework3 R3.5 & R3.11

Homework 3 Problem R*3.5 – Condition for exactness of a N1-ODE
sec13-3

Statement
Prove that
 * {| style="width:100%" border="0"

$$ \bar b(x,y) c(y) y' + a(x)\bar c(x,y) =0 $$     (3.5.1)
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }

can be transformed to an exact form by finding out the integrating factor $$ h(x) $$, only if  $$ k_1(y) = constant$$

Given

 * {| style="width:100%" border="0"

$$ \bar b(x,y) := \int^x b(s)ds + k_1(y) $$      (3.5.2)
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }




 * {| style="width:100%" border="0"

$$ \bar c(x,y) := \int^y c(s)ds + k_2(x) $$      (3.5.3)
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }

Solution
From Eq. (3.5.1), we can see that
 * {| style="width:100%" border="0"

$$ M(x,y) = a(x)\bar c(x,y) $$      (3.5.4)
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }


 * {| style="width:100%" border="0"

$$ N(x,y) = c(y)\bar b(x,y) $$      (3.5.5)
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }

If we assume that $$ h $$ is a function of x only,

From (2)p11.2 we see that the RHS should be a function of x alone. Evaluating $$ N_x $$ and $$ M_y $$,


 * {| style="width:100%" border="0"

$$ N_x = \frac{\partial}{\partial x} \bar b(x,y) c(y) $$      (3.5.6)
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }

From 3.5.2, using the expression for $$ \bar b(x,y) $$


 * {| style="width:100%" border="0"

$$ N_x = c(y) b(x) $$      (3.5.7)
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }


 * {| style="width:100%" border="0"

$$ M_y = \frac{\partial}{\partial y} \bar c(x,y) a(x) $$      (3.5.8)
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }

From 3.5.2, using the expression for $$ \bar c(x,y) $$


 * {| style="width:100%" border="0"

$$ M_y = c(y) a(x) $$      (3.5.9)
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }

Putting these in (2) 11.2,


 * {| style="width:100%" border="0"

$$ \frac{h_x}{h}= -\frac{1}{N} (N_x-M_y) $$      (3.5.10)
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }


 * {| style="width:100%" border="0"

$$ \frac{h_x}{h}= -\frac{1}{\bar b(x,y)c(y)} (c(y)b(x)-a(x)c(y)) $$      (3.5.11)
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }

Cancelling out c(y) from the equation,
 * {| style="width:100%" border="0"

$$ \frac{h_x}{h}= -\frac{1}{\int^x b(s)ds + k_1(y)} (b(x)-a(x)) $$      (3.5.12)
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }

It can be seen that for the right hand side to be a function of only x, $$ k_1(y) $$ should be a constant.

Author and References

 * Solved and Typed by -- Pushkar Mishra

Homework 3 Problem R*3.11 – System of coupled pendulums
sec15-5

Statement
Given a set of coupled equations describing the motion of two pendulums,


 * {| style="width:100%" border="0"

$$ m_1 l^2 \ddot \theta_1 = -ka^2 (\theta_1 - \theta_2)-m_1 g l \theta_1 + u_1 l $$ (3.11.1)
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }


 * {| style="width:100%" border="0"

$$ m_2 l^2 \ddot \theta_2 = -ka^2 (\theta_2 - \theta_1)-m_2 g l \theta_2 + u_2 l $$ (3.11.2)
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }

1. Use ode45 command in MATLAB to integrate the above system by putting it in a matrix form, and integrating for t=[0,7].

2. Use (2) p.15-2, i.e. the solution using IFM for the above equations.

3. Plot $$ \theta_1(t) $$ and $$ \theta_2(t) $$ from Q1 and Q2

Given

 * {| style="width:100%" border="0"

$$ a=0.3 ,\ l=1, \ k=0.2 ,\ m_1g=3, \ m_2g=6 $$      (3.11.3)
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }

Solution
Putting the values of a, k ,m and l in Eq. 3.11.1 and 3.11.2 ,

we get,


 * {| style="width:100%" border="0"

$$ m_1 l^2 \ddot \theta_1 = -ka^2 (\theta_1 - \theta_2)-m_1 g l \theta_1 + u_1 l $$ (3.11.4)
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }


 * {| style="width:100%" border="0"

$$ m_2 l^2 \ddot \theta_2 = -ka^2 (\theta_2 - \theta_1)-m_2 g l \theta_2 + u_2 l $$ (3.11.5)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Dividing Eq. 3.11.1 through out by $$ m_1 l^2 $$ and Eq. 3.11.2 by $$ m_2 l^2 $$
 * {| style="width:100%" border="0"

$$ \ddot \theta_1 = \frac{-ka^2}{m_1l^2}\theta_1 + \frac{-ka^2}{m_1l^2}\theta_2 - \frac{m_2gl}{m_2l^2}\theta_1 + \frac{u_1l}{m_1l^2} $$      (3.11.6)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


 * {| style="width:100%" border="0"

$$ \ddot \theta_2 = \frac{-ka^2}{m_2l^2}\theta_2 + \frac{-ka^2}{m_2l^2}\theta_1 - \frac{m_2gl}{m_2l^2}\theta_2 + \frac{u_2l}{m_2l^2} $$      (3.11.7)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Cancelling and rearranging terms gives,


 * {| style="width:100%" border="0"

$$ \ddot \theta_1 = (\frac{-ka^2}{m_1l^2}-\frac{g}{l}) \theta_1 + \frac{ka^2}{m_1l^2}\theta_2+ \frac{u_1}{m_1l} $$      (3.11.8)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }


 * {| style="width:100%" border="0"

$$ \ddot \theta_2 = \frac{ka^2}{m_1l^2} \theta_1 + (\frac{-ka^2}{m_1l^2}-\frac{g}{l})\theta_2+ \frac{u_2}{m_2l} $$      (3.11.9)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

To put the above equation in the form ,


 * {| style="width:100%" border="0"

$$ \dot x(t) = A(t)x(t)+B(t)u(t) $$      (3.11.10)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

We should define


 * {| style="width:100%" border="0"

$$ x := [\theta_1 \ \dot\theta_1 \ \theta_2 \ \dot\theta_2]^T $$      (3.11.11)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Eq. 3.11. and 3.11. can be rewritten in the matrix form as,
 * {| style="width:100%" border="0"

$$ \begin{Bmatrix}
 * style="width:95%" |
 * style="width:95%" |

\ \dot \theta_1 \\ \ddot \theta_1 \\ \dot \theta_2 \\ \ddot\theta_2

\end{Bmatrix} = \begin{Bmatrix} 0 & 1 & 0 & 0\\ (\frac{-ka^2}{m_1l^2}-\frac{g}{l}) & 0 & \frac{ka^2}{m_1l^2} & 0\\ 0 & 0 & 1 & 0\\ \frac{ka^2}{m_2l^2} & 0 & (\frac{-ka^2}{m_2l^2}-\frac{g}{l}) & 0 \end{Bmatrix} \begin{Bmatrix} \theta_1 \\ \dot \theta_1 \\ \theta_2 \\ \dot \theta_2

\end{Bmatrix} + \begin{Bmatrix} 0 & 0\\ \frac{1}{m_2l} & 0\\ 0 & 0\\ 0 & \frac{1}{m_1l} \end{Bmatrix} \begin{Bmatrix} u_1\\ u_2 \end{Bmatrix}

$$      (3.11.12)
 * <p style="text-align:right">
 * }

Using the values of k,a,m,l and u, the matrix becomes,


 * {| style="width:100%" border="0"

$$ \begin{Bmatrix}
 * style="width:95%" |
 * style="width:95%" |

\ \dot \theta_1 \\ \ddot \theta_1 \\ \dot \theta_2 \\ \ddot\theta_2

\end{Bmatrix} = \begin{Bmatrix} 0 & 1 & 0 & 0\\ -10.06 & 0 & 0.06 & 0\\ 0 & 0 & 1 & 0\\ 0.03 & 0 & -10.03 & 0 \end{Bmatrix} \begin{Bmatrix} \theta_1 \\ \dot \theta_1 \\ \theta_2 \\ \dot \theta_2

\end{Bmatrix} $$      (3.11.13)
 * <p style="text-align:right">
 * }

The Following MATLAB program is written to solve the ODE.


 * {| style="width:100%" border="0"


 * style="width:95%" |
 * function main
 * Xo=[1 -2 -0.5 1];
 * tspan = [0,7];
 * tspan = [0,7];


 * [t,X]= ode45(@Funct, tspan, Xo);


 * end


 * function [xdot]= Funct(t,X)


 * A=[0 1 0 0; -10.06 0 0.06 0;0 0 1 0;0.03 0 -10.03 0];


 * xdot = A*X;
 * end


 * <p style="text-align:right">
 * }

2) Using the solution as discussed in R*3.10 ,


 * {| style="width:100%" border="0"

$$ x(t) = exp (A(t-t_0))x(t_0) $$      (3.11.14)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

since u is 0 for our case.


 * {| style="width:100%" border="0"

$$ \begin{Bmatrix}
 * style="width:95%" |
 * style="width:95%" |

\ \dot \theta_1 \\ \ddot \theta_1 \\ \dot \theta_2 \\ \ddot\theta_2

\end{Bmatrix} = exp (A(t))\begin{Bmatrix}

\ 1 \\ -2 \\ -0.5 \\ 1

\end{Bmatrix} $$      (3.11.15)
 * <p style="text-align:right">
 * }

We wrote the following MATLAB code to find the solution,


 * {| style="width:100%" border="0"


 * style="width:95%" |
 * function main1
 * function main1


 * A=[0 1 0 0; -10.06 0 0.06 0;0 0 1 0;0.03 0 -10.03 0];
 * x0=[1 -2 -0.5 1];
 * dx = zeros(length(t), 4);
 * for i=1 : length(t)
 * dx(i,:) = exp(A*t(i))*x0';
 * end


 * To plot the the different solutions we wrote the following code :
 * hold on
 * figure(3)
 * hold on
 * plot (t,X(:,1))
 * plot (t,dx(:,1), 'r-')
 * hold off
 * grid on
 * legend('Q1','Q2');
 * ylabel('x'); xlabel('t'); title('theta1');


 * figure(4)
 * hold on
 * plot (t,X(:,3))
 * plot (t,dx(:,3), 'r-')
 * hold off
 * grid on
 * legend('Q1','Q2');
 * ylabel('x'); xlabel('t'); title('theta2');
 * end


 * <p style="text-align:right">
 * }

3) The graph of $$\theta_1$$ and $$ \theta_2 $$ as obtained from Q.1 and Q.2 are shown below.



Author and References

 * Solved and Typed by -- Pushkar Mishra