Newton Raphson derivation

Discussion in 'Math' started by kokkie_d, Jan 21, 2011.

  1. kokkie_d

    Thread Starter Active Member

    Jan 12, 2009
    72
    0
    Hi,

    I want to calculate "a" in the following equation using the Newton- Raphson method:

    <br />
M = \phi _2*S*\phi _1<br />
<br />
|eig( M)| - 0.8241 = 0<br />

    With M, phi1 nad phi2 are matrices and S is given by:

    <br />
S = \left[ \begin{matrix} 1 & 0 \\ <br />
\frac{Vi}{x2-x1-aV\omega cos(\omega bT) - C} & 1 \end{matrix} \right]<br />

    Besically all the variables are known apart from "a" thus using newton raphson we can find "a".
    The way I have approached his is by taking the derivative of S wit hrespect to "a", thus
    <br />
S' = \left[ \begin{matrix} 1 & 0 \\ <br />
\frac{Vi*V\omega cos(\omega bT)}{(x2-x1-aV\omega cos(\omega bT) - C)^2} & 1 \end{matrix} \right]<br />

    <br />
M' = \phi _2*S'*\phi _1<br />
|eig( M')|  = 0<br />

    In matlab, I enter all these equations and state that for Newton Raphson:
    <br />
a_{n+1} = a_n - \frac{y(a_n)}{y'(a_n)}<br />

    <br />
y(a_n) = |eig(M)| - 0.8241<br />
    <br />
y'(a_n) = \frac{eig(M')}{|eig(M')|}<br />

    I then loop this a couple of times. Every loop generates a new "a" which is entered and the equations calculated again and again derived.

    But you guessed, it is not working. Now, I was hoping someone could tell me if I had made errors towards the derivation during the various steps or not. Or maybe someone spots an other problem?
     
  2. Georacer

    Moderator

    Nov 25, 2009
    5,142
    1,266
    For starters, the derivative of S is not
    <br />
S' = \left[ \begin{matrix} 1 & 0 \\ <br />
\frac{Vi*V\omega cos(\omega bT)}{(x2-x1-aV\omega cos(\omega bT) - C)^2} & 1 \end{matrix} \right]<br />

    but
    <br />
S' = \left[ \begin{matrix} 0 & 0 \\ <br />
\frac{Vi*V\omega cos(\omega bT)}{(x2-x1-aV\omega cos(\omega bT) - C)^2} & 0 \end{matrix} \right]<br />
    .

    Maybe you should check if this solves the problem first, or I 'm gonna have to do some reading before I answer you with confidence.

    I 'm not sure, but can you write that: M'=\Phi 1 \cdot S' \cdot \Phi 2. That would be the second thing I would check, if the first didn't work.
     
  3. kokkie_d

    Thread Starter Active Member

    Jan 12, 2009
    72
    0
    Thanks, I had just arrived at the same conclusion for the 1's in the matrix do not remain.
    Still, does not work, though.

    The next part:
    M'=\Phi 1 \cdot S' \cdot \Phi 2
    I assumed, this can be treated as any constant times a function remains there after derivation.
    y(x) = 5*a*x --> y'(x) = 5*a
    5*y(x) = a*x --> 5*y'(x) = a --> y'(x) = 5*a

    I am fairly confident it is correct (off course there are smarter people than me out there). What I done to test it is to take the derivative of S and then multiply phi2 with the result and then the outcome of that with phi1. while using matrices of letters.
    I then multiplied phi2 with the original S and that result with phi1 and took the derivative of the final. Both approaches arrived at the same result.

    It is the parts after that I am not confident about.
     
  4. Georacer

    Moderator

    Nov 25, 2009
    5,142
    1,266
    The more I look at your problem, the more I find more inaccuracies and problematic spots.


    • Your function for solving is f=|eig(M(a))|-0.8241=0. However, f'=|eig(M(a))|' isn't equal to |eig(M'(a))|, to my understanding. Could you verify this?
    • Furthermore, the method is x_{k+1}=x_k-\frac{f(a)}{f'(a)}, not x_k- \frac{f(a)}{\frac{eig(M(a))}{|eig(M(a))|}}
    • Finally, you haven't verified that \frac{f(a)}{f'(a)} is smaller than b-a, where [a,b] is the domain in which you try to approach the solution
    Can you post the original question, if this comes from a homework problem?
     
  5. kokkie_d

    Thread Starter Active Member

    Jan 12, 2009
    72
    0
    This is exactly that I am trying to find out.
    According to the chain rule the derivate of the absolute of a fuction is the function divided by it's absolute value.
    <br />
{\frac{eig(M(a))}{|eig(M(a))|}}<br />

    The equation:
    |eig(M(T,0,x(0)))| - 0.8241 = 0
    comes from a paper on stability in buck converters.
    I am replicating it to learn from. This is the last equation I need to solve. I understand what they are trying to do in the paper: The eigenvalues (floquet multipliers) of the monodromy matrix (M) indicate whether a system is stable. If the absolute value of the eigenvalue (I know you have two with one approaching zero and the otherone moving away - becoming more negative and its this last one that is used) is outside the circle with radius 0.8241 then the system is unstable. The value of "a" can be manipulated to bring the system back into stability.
    S stands for the Saltation matrix which describes the point where the switching manifold occurs when the system is non smooth, such as in buck converters.

    In the above analysis I tried to keep things simple but I have also attached a document which describes it more in depth. The calculations that is.

    While analysing all this, I noticed I might have done something wrong with relation to the derivation of the absolute value. My derivation initially would be correct, but since its not a plain x I am trying to derive there might be more to it. In the mean time I am open to suggestions.
     
  6. Georacer

    Moderator

    Nov 25, 2009
    5,142
    1,266
    I 'll insist on that one:
    |f(x)|'=\frac{d|f(x)|}{dx}=\frac{d|f(x)|}{df(x)} \cdot \frac{df(x)}{dx}=\pm f'(x)

    That is the correct way to apply the chain rule, in my opinion.

    As for the rest of the paper, it's a bit complex to me. I 'll need more time to decipher it. If I manage to do it.
     
  7. kokkie_d

    Thread Starter Active Member

    Jan 12, 2009
    72
    0
    I am not clear on what you do in the last step.
    The way I understand it is that an absolute function is composite and thus:
    y =|f(x)|= \sqrt{f(x)^2}
    Let's say: z =f(x)^2
    y = \sqrt{z} = z^{0.5}, y' = 0.5 \cdot z^{-0.5} = \frac{1}{0.5\sqrt{z}}
    But the chain rule says that you also need to multiply with the derivative of z, thus:

    y' = \frac{1}{0.5\sqrt{f(x)^2}} \cdot 2\cdot f(x) = \frac{f(x)}{|f(x)|

    That is the derivative of the absolute part, now the derivative of the f(x), which is eig[M]?

    Is that the way to approach the total derivation?

    btw, I realy appreciate your help. it's great having someone to bounce ideas of.
     
  8. Georacer

    Moderator

    Nov 25, 2009
    5,142
    1,266
    Heh, you tried to do something there, but I think I found another flaw.
    We have:
    <br />
y=|f(x)|=\sqrt{f(x)^2}\\<br />
z=f(x)^2\\<br />
<br />
y'=\sqrt{z}'\\<br />
=0.5 \cdot z^{-0.5} \cdot z'\\<br />
=0.5 \cdot z^{-0.5} \cdot 2 \cdot f(x) \cdot f'(x)\\<br />
=\frac{f(x)}{|f(x)|} \cdot f'(x)\\<br />
    You see, you can't just stop the chain rule whenever you want. You have to get down to x'=1. This is why you have to continue and take the derivative of f(x) too, as it is a function of x.

    Finally, as |f(x)| is either f(x) or -f(x), the ratio \frac{f(x)}{|f(x)|} equals  \pm 1.
     
  9. kokkie_d

    Thread Starter Active Member

    Jan 12, 2009
    72
    0
    Ok, I appologise for the -I thought I new the chainrule while I did not. I had a read and this is my new approach:

    <br />
y=|eig(M)| - 0.8241\\<br />
<br />
    for arguments sake we ignore the value of 0.8241 since it will disappear in the derivation. And the full version of S is provided earlier.
    <br />
y = |eig(M)| = \sqrt{(eig(M))^2}\\<br />
z = (eig(M))^2\\<br />
<br />
y'=\sqrt{z}'\\<br />
  =0.5 \cdot z^{-0.5} \cdot z'\\<br />
  = \frac{1}{2 \cdot sqrt{z}} \cdot z'\\<br />
  = \frac{1}{2 \cdot sqrt{(eig(M))^2}} \cdot ((eig(M))^2)'\\<br />
z = eig(M)\\<br />
y' = \frac{1}{2 \cdot sqrt{(eig(M))^2}} \cdot ((z^2)'\\<br />
y' = \frac{1}{2 \cdot sqrt{(eig(M))^2}} \cdot 2 \cdot z \cdot z'\\<br />
y' = \frac{1}{2 \cdot sqrt{(eig(M))^2}} \cdot 2 \cdot eig(M) \cdot (eig(M))'\\<br />
y' = \frac{eig(M)}{|eig(M)|} \cdot (eig(M))'\\<br />

    I think thus far there are no problems.
    Now, I am looking at how to calculate eigenvalues which is:

    <br />
M = \left[ \begin{matrix} M_{11} & M_{12} \\ M_{21} & M_{22} \end{matrix} \right]\\<br />
<br />
eig(M) = \left| \begin{matrix} M_{11}-\lambda & M_{12} \\ M_{21} & M_{22}-\lambda \end{matrix} \right| = 0 \\<br />
<br />
eig(M) = \lambda^2 - M_{11} \cdot \lambda - M_{22} \cdot \lambda + M_{11} \cdot M_{22} - M_{12} \cdot M_{21} = 0\\<br />
<br />
eig(M)' = 2 \cdot \lambda - M_{11} - M_{22} = 0\\<br />

    M_{11}, M_{22} both contain the matrix S and thus the value of "a" within S we are deriving to.
    The next step is the derivative of M, for which I end up getting the following matrix - I have tried the following in matlab and found it to be correct (The first value indicates whether it is \phi_1 or \phi_2 and the second value the location in the 2x2 matrix.  \left[ \begin{matrix} 1 & 2 \\ 3 & 4 \end{matrix} \right]\\
    <br />
M = \phi _2 \cdot S' \cdot \phi _1\\<br />
S' =  \left[ \begin{matrix} 0 & 0 \\ S' & 0 \end{matrix} \right]\\<br />
M' = \left[ \begin{matrix} \phi _{22} \cdot S' \cdot \phi _{11} & \phi _{22} \cdot S' \cdot \phi _{12} \\ \phi _{24} \cdot S' \cdot \phi _{11} & \phi _{24} \cdot S' \cdot \phi _{12} \end{matrix} \right]\\<br />

    And we have the derivative of S' from previous posts.

    So, getting back to our original outcome:
    <br />
y' = \frac{eig(M)}{|eig(M)|} \cdot (eig(M))'\\<br />
y' = \frac{eig(M)}{|eig(M)|} \cdot (\lambda^2 - M_{11} \cdot \lambda - M_{22} \cdot \lambda + M_{11} \cdot M_{22} - M_{12} \cdot M_{21})'\\<br />
y' = \frac{eig(M)}{|eig(M)|} \cdot (2 \cdot \lambda - M'_{11} - M'_{22})\\<br />
    The derivatives of M_{11}, M_{22} are not complex functions but derivations using product rule and quotient rule thus
    <br />
y' = \frac{eig(M)}{|eig(M)|} \cdot (2 \cdot \lambda - M'_{11} - M'_{22})\\<br />
M'_{11} = \phi _{22} \cdot S' \cdot \phi _{11}\\<br />
M'_{22} = \phi _{24} \cdot S' \cdot \phi _{12}\\<br />
    And S' was already derived through earlier with respect to "a".

    Does this make sense? I am not confident on this work since the derivation is with respect to "a" and thus the eigenvalue calculation should not be derived in the way I have done, I think.

    I am accepting thoughts.

    Cheers
     
  10. Georacer

    Moderator

    Nov 25, 2009
    5,142
    1,266
    My concern is elsewhere. The eighenvalues are values that if put in λ would nullify the matrix's determinant. They are not polynomials.

    An example taken by Matlab:eig \left( \left[ \begin{array}{cc} 1 & 2 \\ 3 & 4 \end{array} \right] \right) = \left[ \begin{array}{c} -0.3723 \\ 5.3723 \end{array} \right] .
    That said, I don't know if saying eig(M) = \lambda^2 - M_{11} \cdot \lambda - M_{22} \cdot \lambda + M_{11} \cdot M_{22} - M_{12} \cdot M_{21} = 0 is correct.

    However, I don't know why the .pdf tries to substract \left[ \begin{array}{c} -0.3723 \\ 5.3723 \end{array} \right] from 0.8241. I can't figure out which norm it is using.
    This is really perplexing...
     
  11. kokkie_d

    Thread Starter Active Member

    Jan 12, 2009
    72
    0
    I agree with the eigenvalue calculation. the way I am starting to read this is: the eigenvalue calculation does not need to be derived since this is a calculation that finds the eigenvalues by applying a fixed method: A - I\lambda = 0\\ Where A is a matrix and I a unity matrix.

    I am not sure what you mean with your last comment (the one from the .pdf); I mentioned in the document the form of the outcome of the eigenvalue calculation in Matlab which is a 2x1 matrix. So if you would want to divide Y by Y' then you are actually dividing one 2x1 matrix by another 2x1 matrix. Well, dividing matrices is a no go, thus a potential solution is to look at only 1 value in the either matrix and divide those, example: Y(2,1) / Y'(2,1). That was all I was trying to expain.

    I am going to try:
    y = |eig(M)| - 0.8241 = 0\\<br />
y' = \pm 1 \cdot eig(M')\\
     
  12. Georacer

    Moderator

    Nov 25, 2009
    5,142
    1,266
    One value at a time eh? Why not... we don't have many other options. Just be sure to substitute the eigenvalue's value, as a function of a. Not the A-Iλ=0 equation.
     
  13. kokkie_d

    Thread Starter Active Member

    Jan 12, 2009
    72
    0
    Ok. the following yields reasonable results
    y = |eig(M)| - 0.8241 = 0\\<br />
y' = eig(M')\\

    But only when I start at a value which is already stable and close to the intended value. So I guess I have to implement the  \pm 1 factor in some why while detecting whether or not I am below stability or not.
    At least I now can continue with the analysis.

    Thanks very much.
     
  14. Georacer

    Moderator

    Nov 25, 2009
    5,142
    1,266
    Start-point stability is a must for the Newton-Raphson method. Maybe you could find mathematically some local minimums, but I think for the purpose of that exercise where the range of the solution is more or less known, you could do with guessing too.

    I admire your persistence on the problem!
     
  15. kokkie_d

    Thread Starter Active Member

    Jan 12, 2009
    72
    0
    Hi Georacer,
    Here is result, I have matched this outcome to the original paper and the graphs are identical. I have used some matlab programming to keep the search in the right area.
    http://twitpic.com/3tj8gf

    Again thanks for all your help. I understand this math a lot better.
     
  16. Georacer

    Moderator

    Nov 25, 2009
    5,142
    1,266
    Could you post some Matlab code on calculating the Newton Raphson iterations? I would be much obliged.

    Congratulations on a mission well done!
     
Loading...