3 views (last 30 days)

I am working on a code that has a system of ODEs, but I have never worked with systems with ode45. In the part of the code that I have included, S is the arclength (which is basically the time step of this problem), th is the angle (theta) of the graph, R is the x coordinate, and Z is the y coordinate.

When I run the program as shown below, I get simply a matrix full of NaN, even when I change the initial R value to 0.0001 or something.

Any help would be appreciated.

Also, what is the output? I only want to graph R and Z, not theta

function yp=program(S,y)

th=y(1);

R=y(2);

Z=y(3);

dthdS=-sin(th)/R+Z-2*H;

dRdS=cos(th);

dZdS=sin(th);

yp=[dthdS; dRdS; dZdS];

end

[S,Y]=ode45(@program, [0, 1], [0, 0, 0])

Star Strider
on 8 Aug 2012

Edited: Star Strider
on 9 Aug 2012

% dThXYdS(1) = Theta, dThXYdS(2) = R(t) = x, dThXYdS(3) = Z(t)

H = 10;

dThXYdS = @(t,ThXY) [-sin(ThXY(1))/ThXY(2) + ThXY(3) - 2*H; cos(ThXY(1)); sin(ThXY(1))];

x0 = [0.1; 0.1; 0];

Tspan = [0:0.01:2]';

[T ThXY] = ode45(dThXYdS, Tspan, x0);

figure(8)

plot(ThXY(:,2), ThXY(:,3))

xlabel('R(S)')

ylabel('Z(S)')

grid

Except for the axis labels in the plot, I used the variable designation from your previous post (and my previous answer) rather than change it to match your current variable designation.

The reason you are getting a matrix of NaNs is that your initial conditions are [0 0 0]. So R is zero when theta is zero and of course sin(theta) will be zero as well. By convention, (0/0) = NaN. An initial NaN in a recursive calculation such as yours creates a matrix of NaNs.

Star Strider
on 9 Aug 2012

I just checked it to be sure it's the same code I created and ran successfully earlier. I ran it again just now and it worked fine for me. (I'm running it on 2012a, but I doubt that makes any difference.) With axis equal it plots something that looks like a cycloid with about 6 cycles. ThXY is a [201 x 3] double array.

Since I can't reproduce your error, I invite others to run it as well to see if they have problems with it.

Star Strider
on 9 Aug 2012

You might be able to do that using the Events property described in odeset, the function that creates an options structure for the ODE solvers. I've not yet had occasion to use that option, so I encourage you to experiment.

Also, using format shortEng or format shortE will solve your numerical resolution problem.

Opportunities for recent engineering grads.

Apply TodayFind the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!
## 0 Comments

Sign in to comment.