How to solve an ode using MatLabs ode23 or ode45

  • Thread starter Nissen, Søren Rune
  • Start date
  • Tags
    Ode Ode45
It is appreciated.In summary, the conversation discusses the use of MatLab's built-in functions ode23 or ode45 to draw a curve of the current i(t) in a circuit with specific values. The conversation also addresses the creation of a derivative function and the process of converting an nth-order ordinary differential equation to a first order system of equations. Eventually, the conversation concludes with a successful solution to the problem and a thank you for the help provided.
  • #1
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.
 
Last edited by a moderator:
Physics news on Phys.org
  • #2
Nissen said:
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
 
  • #3
D H said:
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?
 
  • #4
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?
 
  • #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?
 
  • #6
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.
 
  • #7
D H said:
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
Code:
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:
[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.
 

1. How do I define my ODE function in MatLab?

In order to solve an ODE using MatLab's ode23 or ode45 functions, you must first define your ODE function using the appropriate syntax. This function should take in the independent variable (usually time) and the dependent variables as inputs, and return the derivatives of the dependent variables. It is important to make sure that the function is vectorized, meaning it can handle arrays of inputs.

2. What is the difference between ode23 and ode45?

Ode23 and ode45 are both numerical integrators used to solve ODEs in MatLab. The main difference between the two is the order of the method used to approximate the solution. Ode23 uses a second-order method, while ode45 uses a fourth-order method. This means that ode45 is generally more accurate, but may also take longer to compute.

3. How do I specify the initial conditions for my ODE?

In order to solve an ODE using MatLab, you must specify the initial conditions for your dependent variables. This is typically done by creating a vector containing the initial values for each variable, and passing it as an argument to the ode23 or ode45 function. It is important to make sure that the order of the variables in the vector matches the order in which they are defined in your ODE function.

4. Can I plot the solution to my ODE using MatLab?

Yes, MatLab has built-in functions for plotting the solution to your ODE. After using ode23 or ode45 to solve the ODE, you can use the "plot" or "plot3" function to visualize the results. It is important to note that the solution will only be plotted for the range of independent variable values that were specified in the initial call to the ode function.

5. What should I do if my ODE does not converge using ode23 or ode45?

If your ODE does not converge using ode23 or ode45, there are several things you can try to improve the solution. First, you can try changing the options for the ode function, such as the error tolerances or the integration time step. Additionally, you can try using a different integration method, such as ode15s or ode23s. It may also be helpful to double-check your ODE function and initial conditions to make sure they are correct.

Similar threads

  • Calculus and Beyond Homework Help
Replies
4
Views
437
  • Engineering and Comp Sci Homework Help
Replies
1
Views
817
  • Engineering and Comp Sci Homework Help
Replies
2
Views
764
  • Introductory Physics Homework Help
Replies
5
Views
3K
Replies
3
Views
1K
  • MATLAB, Maple, Mathematica, LaTeX
Replies
5
Views
2K
  • MATLAB, Maple, Mathematica, LaTeX
Replies
6
Views
2K
  • Calculus and Beyond Homework Help
Replies
1
Views
4K
  • MATLAB, Maple, Mathematica, LaTeX
Replies
1
Views
1K
  • MATLAB, Maple, Mathematica, LaTeX
Replies
8
Views
1K
Back
Top