System of arc-length defined ODEs with ode45

Question

Michael 2012-8-8

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/45598-system-of-arc-length-defined-odes-with-ode45

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])

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Star Strider 2012-8-8

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/45598-system-of-arc-length-defined-odes-with-ode45#answer_55774

编辑：Star Strider 2012-8-9

在 MATLAB Online 中打开

As with your 'How can I graph a "system" of ODEs?' post, I suggest:

        % 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.

4 个评论
显示 2更早的评论隐藏 2更早的评论

Michael 2012-8-9

编辑：Michael 2012-8-9

I guess my error was rounding and scale, because the graph makes a shape, but the tables show otherwise.

This is looking to be similar to the shape I want. Is there a way to limit it so that the solver stops when it starts to go back up (past a half-cycle). Like saying "stop when z(k)<z(k+1)", or do I have to just limit the Tspan?

Star Strider 2012-8-9

编辑：Star Strider 2012-8-9

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.

请先登录，再进行评论。