Ode 45 Numerical Asnwers

Ok, well this may seem a bit stupid but I have no idea how to find a numerical answer from the ode45 function.

Ok I am trying to solve d2x/dt2 = x(1+Asin(wt)

With initial conditions x(0)=0.1 and dx(0)/dt= 0

Where a and w are constants to be given.

This is what I have for a mfile

function xprime=third(t,x);

A=100;
w=20;

xprime(1)=x(2)
xprime(2)=x(1)*(1+A*sin(w*t)

xprime=xprime(:);

This is fine and by using

ts=[1,4*pi];
  x0=[0.1 0];
  [t,x]=ode45(@third,ts,x0);

plot(x(:,1),x(:,2))

I can get a graph which I assume is correct, but I need to find the numerical solution of when A=100 and w=20 (I also need for when w = 40, 60, 80, 100) and I am not sure where I can get them from.

I am sure it is obvious and would appreciated an answer. Also could someone please check my code as I have no idea if it is right.

Thanks.

[Jiro Doke]

One comment about your graph: You say "I assume [it] is correct". Does it look correct? You're plotting dx/dt vs x. Just want to make sure you weren't looking for time plots.

Answer 1

[Jan]

If you need x(0)=0, dx(0)/dt=0, the initial vector must be x0=[0, 0], not [0.1, 0].

If you want to modify the parameter of the function, use the parameter input of ODE45 as described in the help text.

function xprime = third(t, x, w)
A = 100;
xprime = zeros(2, 1);  % Pre-allocate
xprime(1) = x(2);
xprime(2) = x(1)*(1+A*sin(w*t));

The difference between setting yprime to ZEROS at first or to reshape it after the creation is tiny. For larger arrays the difference is dramatic!

Now your main function:

w = 20; % 60, 80, 100
ts = [1, 4*pi];
x0 = [0, 0];

% EDITED: New syntax for parameters as suggested by Jiro:
[t,x]=ode45(@(T,X) third(T,X,w), ts, x0);
plot(t, x(:,1), t, x(:,2));

Does this help?

[Jiro Doke]

@Jan: Did you mean to use an anonymous function for your call to ode45?

w = 40; % 60, 80, 100

[t,x]=ode45(@(T,X) third(T,X,w),ts,x0);

[Jan]

@Jiro: No, I meant using the parameter w as "p1" according to "doc ode45". Therefore I set the options to [] in the ode45 call.

[Jiro Doke]

@Jan: Which version of MATLAB are you using? In the recent releases, the recommended way of passing additional parameters is to use an anonymous function or a nested function (http://www.mathworks.com/help/matlab/math/bsgprpq-5.html). The documentation doesn't have the old syntax anymore:

http://www.mathworks.com/help/matlab/ref/ode45.html

But the old syntax still works.

[Jan]

@Jiro: Doh. I've written my simulations in Matlab 6.5 and compared the source code of ODE45 when I updated to 2009a - but I forgot to read the help text in detail. Thanks, Jiro!

I definitely prefer the old parmeter style, because I call compiled C-functions usually. In addition the external numerical differentiation looks more convenient when I use the parameters as arguments of the integrator. Is using an anonymous function as efficient as the old style?

[George Harrison]

Hmmm ok I suppose this makes sense, so how do I go about getting the numerical answer then?

[Jan]

What is "the numerical answer"? Do you mean the the last value of x? Then use x(end,1) and x(end, 2).

[Jiro Doke]

What do you mean by "numerical answer"?

"x" that you get from ODE45 is the numerical answer. It's actually vectors (2-column matrix) of answers with the time stamp "t". First column of "x" is the position trace and the second column is the velocity trace.

[Jiro Doke]

@Jan: I asked around about your question about the differences. I can't speak to the "efficiency" question, but I've summarized the response in this new question:

http://www.mathworks.com/matlabcentral/answers/1971-when-using-ode45-or-similar-functions-what-is-the-benefit-of-using-anonymous-functions-over-passi

[George Harrison]

I am asked to "find the numerical solution" for the cases where w = 40 60 and 80. So I am not really sure, I assume it will be the first value of X or possible the last. Either way I would still like to know where to find these.

Thanks for the help so far though :)

[Jiro Doke]

Not sure what more we can give you. "t" and "x" _are_ the numerical solution.

Answer 2

[Jan]

If you integrate an ODE, the solution is the trajectory. Therefore the complete array x and t is the "numerical solution".