How to solve an ode using MatLabs ode23 or ode45

[Nissen, Søren Rune]

Homework Statement

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

Homework Equations

That would be the one above.

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.

 

[D H]

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);

Your derivative function is incorrect. It is not a re-expression of

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

 

[Nissen, Søren Rune]

Your derivative function is incorrect. It is not a re-expression of

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?

 

[D H]

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?

 

[Nissen, Søren Rune]

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?

 

[D H]

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 n^th-order ordinary differential equation

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

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

\aligned<br /> y_1 &amp;\equiv y&#039; \\<br /> y_2 &amp;\equiv y&#039;&#039; \\<br /> \cdots \\<br /> y_{n-1} &amp;\equiv \frac{d^{n-1}y}{dx^{n-1}}<br /> \endaligned

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

\aligned<br /> y&#039; &amp;= y_1 \\<br /> y_1&#039; &amp;= y_2 \\<br /> \cdots \\<br /> y_{n-2}&#039; &amp;= y_{n-1} \\<br /> y_{n-1}&#039; &amp;= f(x,y,y_1,\cdots,y_{n-1})<br /> \endaligned

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

The error in your function current is a simple algebraic mistake.

 

[Nissen, Søren Rune]

The error in your function current is a simple algebraic mistake.

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

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

[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.