ode45 problem with plotting

Hi all,

I am trying to run this code, aparently it has no problem. But it should give me some plots which it does not. I would be happy if you could help me out.

Best,

Argu

[madhan ravi]

Why do you use delete([g1 g2]) ?

[Argumanh]

if you dont do that in viewmotion you are gonna see a trace of the pendulum then.

Answer 1

[KSSV]

You should change view to 1. And there is small change in the name of variable..hope the below one works:

function [t,th,r,J,om]=lab1(eps,view)

% user input:
% (1) eps: this is the variable for the pendulum length: r(t)=r0*(1+eps*t); see lab sheet
% (2) view: if user enters view > 0, a window will pop up on screen to show the
% pendulum motion after the simulation is finished. It is recommended to
% set view=0 since the processing time required to show graphics on screen.
%
% program output:
% (3) t: the time array of 500,000 elements, 0, 0.001, 0.002, ... , 500
% (4) th: a 500000x2 array where each row of this array holds the record
% of the pendulum angle, theta, (d/dt)theta at time t. In particular,
% th(:,1) contains 500,000 elements of theta value at the time instants
% defined in the t array; see (3), and th(:,2) contains the 500,000 elements
% of the (d/dt)theta at the time instants defined by the t array.
% (5) r: the pendulum length as a function of time defined in t array.
% (6) J: the approximation of the total energy of the pendulum mass; see lab sheet
% (7) om: the instantaneous frequency of oscillation of the pendulum; see lab sheet

% assume pendulum mass m=1

r0=1;           % pendulum length at t = 0
dt=0.001;       % integration time step
tspan=0:dt:500; % all time instants for integration

th0=[0.2;0];    % initial condition: theta(0)=0.2, (d/dt)theta(0)=0



opt=odeset('Reltol',1e-9,'AbsTol',1e-9);
[t,th]=ode45(@pend,tspan,th0,opt,r0,eps);

r = r0*(1+eps*t); % pendulum length over time

J=0.5*(r.^2.*th(:,2).^2)+0.5*9.81*r.*th(:,1).^2;
om=sqrt(9.81./r);

if view
    viewmotion(t,th,r);
end

    function dthdt=pend(t,th,r0,eps)
        r=r0*(1+eps*t); % pendulum length at time t
        rdot=r0*eps;    % d/dt pendulum length at time t.
        dthdt(1,1)=0;   % ... enter statement for (d/dt) theta
        dthdt(2,1)=(9.81*sin(0.2)-2*rdot*0)/r;   % ... enter statement for (d^2/dt^2) theta
    end

    function viewmotion(t,x,r)
        clf;
        plot([-1,1],[0,0],'color',[1 1 1]*.5,'linewidth',5);
        hold on;
        plot([0,.5],[0,.2],'color',[49 79 79]./256,'linewidth',3);
        axis([-1.5 1.5 -max(r) 1]);
        npt=length(t);
        for i=1:100:npt
            xi=r(i)*sin(x(i,1));
            yi=-r(i)*cos(x(i,1));
            g1=plot([0,xi],[0,yi],'color',[49 79 79]./256,'linewidth',3);
            g2=plot(xi,yi,'r.','markersize',75);
            pause(1e-2);
            delete([g1 g2]);
        end
        
    end

end

[Argumanh]

I would appreciate if you also could attache the file instead of writing it here. I will try the new code.

Thanks,

Argu

[KSSV]

There is a copy option......simply copy and paste......I should go for more mouse clicks to attach....:)

[Argumanh]

I know, but I would highly appreciate if you edit your answer like that. The thing is, the code is working. But I cannot see stored data as explained in the output part of the code(the statement at the first of the code)

best,

Argu

[KSSV]

Edited...I guess this what you want.

[Argumanh]

now it says

Not enough input arguments.

Error in lab111 (line 24)

[t,th]=ode45(@pend,tspan,th0,opt,r0,eps);

Also I would appreciate if you could attach the file. Thanks

[Argumanh]

I want to have that ''th'' array. I do not ge it, it just gives "ans"

[KSSV]

eps = 0.1 ;
view = 0 ;
[t,th,r,J,om]=lab1(eps,view) ;

The above is the way to use the function.

[Argumanh]

then how can I get the "th" array?

[KSSV]

Already you have th array in hand.......see the outputs of the function.

[Argumanh]

I just get "ans" in the workspace it soed not give me the "th" array. I got confused

[Argumanh]

Dear KSSV,

I still do not understand how to get the "th" array. as I said after running the code on the workspace I can see "ans" there is no other arrays on it. How can I get that?

Best