Problem with ode45 code ?

Hi, i'm really new to Matlab. I am using ode45 to solve this differential equation:

y1prime = y2

y2prime = -y1 - (1/t)*(y2)

I am plotting it over the closed interval 1, 50.

It has the intial conditions:

y1(0)=0
y2(0)=0
y2prime(0)=0.

However i keep getting error messages about my m-script. I am using a code from quite an old book, so i'm not sure if that has anything to do with it. I would appreciate if anybody could have a quick look to see where i've gone wrong. Thank you

function yprime = pend(t,y)

yprime = [y(2); -y(1)-((y(2))/t)];

tspan = [1 50];

yazero = [0;0]; ybzero = [0;0]; yczero = [0;0];

[ta, ya] = ode45(@pend, tspan, yazero);

[tb, yb] = ode45(@pend, tspan, ybzero);

[tc, yc] = ode45(@pend, tspan, yczero);

[y1,y2] = meshgrid(-5:.5:5,-3:.5:3);

Dy1Dt = y2; Dy2Dt = -y(1)-((y(2))/t);

quiver(y1,y2,Dy1Dt,Dy2Dt)

hold on

plot(ya(:,1),ya(:,2),yb(:,1),yb(:,2),yc(:,1),yc(:,2))

axis equal, axis([-5 5 -3 3])

xlabel y_1(t), ylabel y_2(t), hold off

end

[John D'Errico]

I seriously doubt that ANY old book has MATLAB code in it where you create a function that calls ODE45, inside passing ODE the name of the function that you just created, to create a recursive nightmare of a loop. So you have seriously hacked around with whatever code you did find.

As well, you state:

"I am plotting it over the closed interval 1,50.

It has the intial conditions: y1(0)=0 y2(0)=0 y2prime(0)=0."

But then specify initial conditions at 0, not 1.

I also have no idea why you think you need to solve it three times, all with the same initial conditions. Would you get different answers each time? No.

Finally, you have only two variables, not three. So specifying three initial conditions will be difficult. Anyway, y2prime appears to be trivially zero IF we assume that you really mean to start at 1 for t, NOT 0.

[Jenny Briggs]

I seriously didn't change or hack around much with the code. I don't know enough to change that much. The orginal code in the book doesn't work on my computer either.

[Jan]

@paige: Whenever you mention, that some code "does not work", post the error message or explain the difference between teh results and your expectations. Solving problems is easier than guessing them.

[Brian Russell]

I see your problem. The old book you refer to is "MATLB Guide" by Higham and Higham, which I think is in general a brilliant book. However, they often leave out statements in their code which they assume their reader will know to put in. I also struggled with this example until I realized that their function is very short. In fact, it is only two lines:

function yprime = pend(t,y)

yprime = [y(2); -y(1)-((y(2))/t)];

There should be an end after these two lines. This short function should be saved in a file called pend.m. The rest of your code should be put in another program, say pend_run.m. This will fix your problem.

[Brian Russell]

I meant "MATLAB Guide", of course!

Answer 1

[John D'Errico]

So, just suppose I decide to try to solve your problem, assuming that you really intended to solve it over the interval [1,50]. That would take no more complicated code than this:

yp = @(y,t) [y(2);-y(1) - y(2)./t);
[tout,yout] = ode45(yp,[1,50],[0 0]);

The result would be rather boring though, since you have postulated a differential equation system of the form:

y1' = 0
y2' = 0

Note that since y1(1) = y2(1) = 0 at the beginning, they will STAY at zero forever. And ever. And ever.

Kind of boring.