Solving Second Order Differential Equations in Matlab

[smithforthphd]

thanks! I am a professor here at Oregon State. I am going to take this to the classroom!

[danielalonso23]

Thanks you very much dude¡¡¡¡¡¡, Greetings from Chile

[alpha001ful]

How can define " a function that defines the functions derived on previous slide (the first deivatives of Z and Y)". i think it is not shown here and that is my problem.

[jakeblanchard]

@alpha001ful

It is easiest if you send me a script I can debug. Send it to blanchard "at" engr.wisc.edu

Alternatively, download a sample script from

blanchard "dot" ep "dot" wisc "dot" edu/PublicMatlab/#IVP

[alpha001ful]

I tried to do the same simulation for a simple pendulum. Instead of y there , i had angle or variable called theta. Whenever i tried to run, Matlab always said "??? Input argument "theta" is undefined."........i have done exactly the same as here...

[tomtomtramp]

I've been using your videos, thanks!. But, I think I'm writing the wrong code. I am developing a system to analysis the rocking motion of a block subjected to an earthquake excitation. The state vector is y(t)={theta(t) , theta'(t)}. and the Time-derivative vector f(t) is y'(t)= {theta'(t) , p^2(sin(alpha*sgn[theta(t)]-th­eta(t))+u''/g*cos(alpha*sgn[th­eta(t)]-theta(t))]}. I have the values for p, alpha, u'' and g. But I need theta(t), and theta'(t). Can anyone help, I really need it!

[jakeblanchard]

@tomtomtramp

You can email me at blanchard "at" engr.wisc.edu

[dequinox]

Thanks for this, I couldn't figure out why my code wasn't working...everything looked good. Then I realized I forgot to define g. Dang I want those two hours of my life back! Lol

[pikachu832]

I have a question, if it is not obvious that you are given a 2nd order diffrential equation, how do you go about making it into one? Also, the last matrix at about 5.25; do they have to be [dy/dt; d2y/dt2]??? Many thanks.

[jakeblanchard]

@pikachu832

You need to end up with two first order equations. A 2nd order equation can easily be converted into a pair of first order equations. I'm not sure what else you might run into to make you wonder what you are working with. I'd have to see an example.

The term [dy/dt;d2y/dt2] is actually the first derivatives of the two first order ODE's. If you start with a second order equation, then this notation suggests one possible way to convert to first order equations, but there are others.

[pikachu832]

@jakeblanchard I think you have cleared a lot of things for me. =) But just to confirm... how does matlab know that y(2) = dy/dt? Why is it not y(1)? I suspect it has something to do with the initial values you specified later; if y(2) = dy/dt, then y(2) =z. The 2nd variable in your initial values so happens to be z0 (zinitial). So can I ask if there is a direct link between the inital values implemented and the vector y that MATLAB reads in the previous function file? Thanks again. =)

[jakeblanchard]

@pikachu832

I start with the first question. Matlab doesn't really know y(2)=dy/dt. We decide to have dy/dt=z and we let z=y(2). So our first equation is dy/dt=z. Then, dz/dt=d2y/dt2 and we can find an equation for dz/dt in terms of y and z from the original second order equation. This gives us our second first order ODE. Then we tell Matlab what those two equations are. That's all it knows. This is all completely independent of the initial values.

[jakeblanchard]

@pikachu832

Now for the second question. The differential equations are independent of the initial conditions. When Matlab solves the system, it starts with the initial conditions and then steps its way through time using the differential equations.

[davidwooding]

please please please help me as everyone in my class it stuck on this question as it makes no sense to us.

Part (3) Integrating p'' + 0.1536*p' + k^2*p = 0 with intial conditions p(0)=0, p'(0)=1, plot p(L,k) as a function of k.

[davidwooding]

@davidwooding In this case, you are given the initial conditions (p(0)=0, p'(0)=1), but you must guess a value for parameter (k) by trying a sample of values for k.

For each k-value guess, you can calculate a solution p(x,k) and match this to see how good this is against the end condition (x=L) for the original BVP, p(L)=0. Hence, the problem is all about finding the correct k's, so that in the mode solutions p(L,k)=0

[jakeblanchard]

@davidwooding

I'm not sure what your question is. What you describes seems workable. It looks to me like it is using the shooting method for an eigenvalue problem.

If you guess k (k=1, for example), then you should be able to solve the problem using the approach described in this video. Then you repeat for many values of k and look for the value of k that gives you p=0 at x=L.

[davidwooding]

@jakeblanchard the next part of the question is using the shooting method find resonance wavenumbers k with relative accuract 10^-4 find resonance frequencies f. When you say look for the value of k which gives p=0 at x=L how do i do this? I have everything working from your method :) but this part is really confusing me. Thank you so much for your help.

[jakeblanchard]

@davidwooding

As you said above, can you just plot p(L) as a function of k and see if it ever crosses 0?

[davidwooding]

@jakeblanchard i am really sorry but this is where i am confused. I have plotted plot(x,p(:,1)) for what ever value of k i chose the graph will always cross the x axis where y=0 Im really grateful for your help, i appreciate this most likely doesnt come across well and might be confusing to see what i have done without seeing my mfile

[jakeblanchard]

@davidwooding

You don't want to plot p vs x. You want to just calculate the last value of p. Something like p(end,1). Then you pick a new value of k and rerun the whole thing to, again, find p(end,1). Do this many times, collect all the values for p(end,1), and then plot those vs. the values of k used to generate them.

[davidwooding]

@jakeblanchard thank you so much for all your help it was really appreciated :)

[RightRocks]

I'm sorry, one more question,

What if there is a term involving the independent variable multiplying one of the derivatives?

Say, d2y+x*dy=0? Would I just multiply p(2) by x?

[jakeblanchard]

@RightRocks

Yes, if you have a term like x*(dy/dx) it's just x*y(2)

[RightRocks]

Ah, I get it now. Thank you very much.

[RightRocks]

At 5.28 shouldn't the second y be y(1) rather than y(2) if it represents velocity?

[jakeblanchard]

@RightRocks

No. y(1) is the height (y). y(2) is the first derivative of the height, that is, the velocity. So the second y should be y(2) because the drag term depends on v^2.

[RightRocks]

@jakeblanchard

Ah I see. So with the equation d2y/dx2+m*dy/dx+c*y, I would write the function as:

rk=[y(2), -m*y(2)-c*y(1)];

Is that right?

[jakeblanchard]

@RightRocks

Yes, that is correct.

[baderb52]

THNK you thhhhaaaaaaaaaaaaaaaaaank you you saved my life

[wildreams]

have a question.. shouldn't it be 4/15 instead of 4/14 in your ode function?

[jakeblanchard]

@wildreams

Yes, you are correct. It should be 4/15

[golfer1059]

great vid thanks a bunch!

[bancuk]

@jakeblanchard

what if the coefficient of the ODE is a variable:

d²T/dx² -(m+n sin(kx)) T = q + p sin (kx)

m,n,q, and p are constant

how to solve this?

[jakeblanchard]

@bancuk

This is pretty straightforward. The function that defines the equation will look something like this:

function rk=f(x,y)

y will be a vector of T and dT/dx. As you can see, the current value of x is there in the argument list, so you can use it to define the equations. So you would have something like this

rk=[y(2), (m+n*sin(k*x))*y(1)+q+p*sin(k*­x)];

[bancuk]

@jakeblanchard

what does y(2) and y(1) means? I have never used matlab to solve differential equation? How can y has two meanings, as T and dT/dx?

Thanks.

[jakeblanchard]

@bancuk

When you solve a second order differential equation, you are actually solving two first order equations. So there are two results, one for y and the second for dy/dx. So Matlab saves both in a vector. Typically, y(1) would represent T and y(2) would represent dT/dx.

By the way, you may be solving a boundary value problem, rather than an initial value problem. These require a somewhat different approach.

[topito2]

thank you very much, nice videos!!!!! very informative

[vnck25]

Was Matlab invented by one Dr. Matthew Labvantos? Thus the name MatLab.

[jakeblanchard]

@vnck25

:)

Actually, I think it came from Matrix Laboratory

[vnck25]

@jakeblanchard

Probably thats where this Dr. Labvantos works or worked right?

[bangthatdrumb]

thank u

[pawanmandawat]

Hey Just a quick question, how can we call this 'f.m' file using

[t,y]=ode45(@secode,[0 30],[600 0]); in another *.m file because here what we are doing is calling it in Command Video.

Because when i am trying it in a *.m file it is giving me the following error...

??? Undefined function or method 'falling' for input arguments of type 'char'. where falling is my file name.

Thanks a lot.....

[jakeblanchard]

@pawanmandawat

I'm not sure how to answer this because you refer to f.m, secode.m, and maybe falling.m. Let me try something else. Look in the next posting. I will paste in a script that will work.

[jakeblanchard]

@jakeblanchard

OK. That isn't working. YouTube must think I'm posting spam, or a URL, or something. There is a script called ode2ndOrder.m at my web site blanchard "dot" ep "dot" wisc "dot" edu/PublicMatlab

Download this script and run it and it should give you a working solution.

[pawanmandawat]

Yeah sorry about that,i didn't mentioned that, yeah the file 'f.m' here is 'secode.m' and what i basically want is to call this 'f.m' or 'secode.m' function in another '*.m' file which is 'falling.m' not in command wiindow.

Thanks very much for reply.

[pawanmandawat]

Hey i got the way to do that... thanks very much for showing the way.......

thanks a ton......

waiting for your other videos..........

[BLkMan2TaLk2]

Helpful stuffs ... I appreciate this, man.

[DotComMillie]

Thanks very much, this was a big help. =)

[10linard]

how can use the command diff in matlab in this function...x^3+x^2-1...

[jakeblanchard]

If you have the symbolic toolbox, you can do this

x = sym('x');

diff(x^3+x^2-1)

The result is

ans =

3*x^2 + 2*x

[kidwhocracked]

Thank you!

[parthsna]

Legend

[saadahmedmahdi]

Good explain (Y)

Thank you

[WorldEditna]

hi again. im really lost now with my differential equation (d2i/dt2)=(sin(w*t)/(L*C*Rs))-­(di/dt)/(L*Rs)-i/(L*C).

i need to solve it with Eulers method and i dont understand what should i enter to matlab to solve.

[jakeblanchard]

First you break the second order equation into two first order equations. Then you put together an Euler method for the system. Wikipedia has some info on this.

[kksenduran]

jhj

[samarpan88]

I tried to follow the instructions given here but i m getting the following error-

Strings passed to EVAL cannot contain function declarations.

Is this bcoz. i m using a old version(3.6.5)?

Plz help

[jakeblanchard]

It could be a version problem. Old versions of Matlab didn't use @f. Try using 'f' instead of @f. That is, try single quotes around the f, like this

[t,y]=ode45('f',[0 15], [600,0]);

[zoro19840]

I treied it, it doesnt work ,

how you run this ?

you create 2  .m files ? where you saved them what ar ethe names ?

coul you please answer

[jakeblanchard]

Put these commands into a file called f.m

function rk=f(t,y)

rk=[y(2); -9.81+4/14*y(2).^2/80];

Then type these into the command window

[t,y]=ode45(@f,[0 15], [600,0]);

plot(t,y(:,1))

[zoro19840]

Thank you very much !!

Now it works,

I need another favor please:

how to solve : m(d^2/dt^2)+c(dx/dt)+kx=u(t);

with initial conditions: x(0)=1, (dx/dt)(0)=2; and u(t) is input = sin(2*pi*t)

[jakeblanchard]

Take a look at this site:

blanchard "dot" ep "dot" wisc "dot" edu/PublicMatlab/index.html#IV­P

[WorldEditna]

or juste show how to enter your ode to matlab

[WorldEditna]

Hi! can u help me to solve a second differantial equation based on LC circuit, using matlab, here it is : d^2(I)/dt^2={-1/sqrt(L*C)}*I. thanks

[jakeblanchard]

First, you need a function that defines the equation. Something like this:

function d=rklc(t,y)

i=y(1); z=y(2);

L=1; C=1;

d=[z; -1/sqrt(L*C)*i];

[jakeblanchard]

Then, to get a solution, you do something like this:

init=[1 0];

time=[0 10];

[t,y]=ode45(@rklc, time, init);

plot(t,y(:,1))

[WorldEditna]

thanks for ur help

[diegoramos27]

I hate matlab it's so complex, good tutorial tho.

[Sean6196]

lol well for some people matlab is like a lifeline....for example aerospace engineers :P

[clustro]

nub

[JoeyA73]

One word. simulink.

[zoro19840]

Simulink for kids !!!

[asakuyoshida]

excellent video!

It's very useful.

thanks!

[AntiGravity911]

This video is amazing ! Even though youtube has many dumb videos, it is an excellent resource for education materials

[pijan52]

Good video: short and to the point.

[paul5114]

Great video man! It really helped me with part of my coursework, anyway think there could be an error in your code...

In the calculations, you put 4/15 as the constant but in the code, you put 4/14...Is that just me missing something?

[jakeblanchard]

Good catch. You are correct. They should both be the same.