1. Not finding help here? Sign up for a free 30min tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

How to solve an ode using MatLabs ode23 or ode45

  1. Nov 20, 2008 #1
    1. The problem statement, all variables and given/known data

    Use one of MatLabs built in functions ode23 or ode45 to draw a curve of the current i(t) when

    0.001i''(t) + 200i'(t) + 1/(20*10^-9)i(t)=0

    and

    i(0)=0
    i'(0)=10000

    2. Relevant equations

    That would be the one above.

    3. The attempt at a solution

    I know from my work doing this by hand that the solution is

    i(t)=(1/20)*e^(-100000t)*sin(200000t)

    so I limit my t-span to [0,0.001] - any more and the curve would be invisible.

    Doing this in MatLab:

    The file current.m has been created, and currently contains the following:

    function di=current(t,i)
    di=zeros(2,1);
    di(1)=i(2);
    di(2)=-2*(10^5)*i(2)+(2*10^-7)*i(1);

    The command given in the commandline window to create the I(T) values is

    [T,I]=ode45(@current,[0,0.001],[0,100000])

    followed by plot(T,I)

    This gives me an curve that describes an exponential fall from 10^5 to 10^4 over the span of 10^-4 seconds.

    This is completely at odds with what I'd expect, looking at the curve I calculated by hand (and completely at odds with an underdampened RLC circuit, which is what the original equation describes - 1mH, 200 ohm and 20 nanofahrad)

    EDIT: I should probably point out that I is returned in complex numbers - this is not wholly unexpected, and I suspect that if I knew how to turn the complex I into a sinus rythm, I might get the curve I expected in the first place. I should probably also point out that it's one in the AM where I'm at, so I'm heading to bed. I'll check this thread again tomorrow to look for suggestions or add clarifying information if any is needed.
     
    Last edited by a moderator: Nov 20, 2008
  2. jcsd
  3. Nov 20, 2008 #2

    D H

    User Avatar
    Staff Emeritus
    Science Advisor

    Your derivative function is incorrect. It is not a re-expression of
     
  4. Nov 20, 2008 #3
    Could you tell me whether the error is in line 1, line 2 or in both lines?

    Yes I know I said I'd get to bed - I just can't sleep when this is rummaging around in the back of my brain, apparently.

    EDIT: I'm not surprised the error is in that part - I'm not even sure what "di(1)=i(2)" means in this context. I get that it replaces the first zero in the di matrix with the second number in the i matrix - but what is the i matrix?
     
  5. Nov 20, 2008 #4

    D H

    User Avatar
    Staff Emeritus
    Science Advisor

    I'm not going to tell you what's wrong; that is a problem for you to solve.

    What you need to do is to re-express that second-order equation in one variable to a first-order equation in two variables, typically by defining an auxiliary variable as the first derivative of the variable of interest. You have done just this, using the i array is to represent i and di/dt. Your derivative function, in this case current needs to return the derivatives of the components of your state vector.

    So, what are the time derivatives of i and di/dt, expressed in terms of i and di/dt?
     
  6. Nov 21, 2008 #5
    I think I've found where I'm doing wrong: I need to rewrite my original equation into Cauchy form - but I've received no training in how to do this. Is there a good tutorial on how to do this somewhere on the net?
     
  7. Nov 21, 2008 #6

    D H

    User Avatar
    Staff Emeritus
    Science Advisor

    All you need to do is to reduce the problem to a first-order system of equations. I don't know if this is what you mean by "Cauchy form"; that term doesn't generate many hits on an internet search.

    The procedure for converting an nth-order ordinary differential equation

    [tex]\frac{d^ny}{dx^n} = f(x,y,y',...y^{(n-1)})[/tex]

    to a first order system of equations is to introduce n-1 variables that represent the first n-1 derivatives of the original variable:

    [tex]\aligned
    y_1 &\equiv y' \\
    y_2 &\equiv y'' \\
    \cdots \\
    y_{n-1} &\equiv \frac{d^{n-1}y}{dx^{n-1}}
    \endaligned[/tex]

    The derivatives of the n variables y, y1, ... yn-1 expressed in terms of these variables are thus

    [tex]\aligned
    y' &= y_1 \\
    y_1' &= y_2 \\
    \cdots \\
    y_{n-2}' &= y_{n-1} \\
    y_{n-1}' &= f(x,y,y_1,\cdots,y_{n-1})
    \endaligned[/tex]

    For a second-order equation you need to introduce but one variable.

    The error in your function current is a simple algebraic mistake.
     
  8. Nov 21, 2008 #7
    Yes. Yes it is. One of the people in my group had found it when we returned to class today.

    So we got it to work (finally)

    Here's the function:

    current.m
    Code (Text):
    function di=current(t,i)

    %0.001*i''+200*i'+(1/2*(1/(20E-9)))*i=0

    di=zeros(2,1);

    di(1)=i(2);

    di(2)=-2*(10^5)*i(2)-(5*10^10)*i(1);
    current2.m
    Code (Text):
    [T,I]=ode23(@current,[0,0.0001],[0,10000]);

    plot(T,I(:,1))
    Paints a nice curve looking exactly as expected. Finally.

    Thank you for your help.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?



Similar Discussions: How to solve an ode using MatLabs ode23 or ode45
Loading...