Using ode45 to solve a Second Order Non-linear ODE

I'm trying to solve a part of the following equation :

The part i'm trying to solve is or

I've rewritten itt in the form of two first order equations, like so:

In vector form,

And my code is as follows:

tspan = 0:0.0033:100;
a=10*(pi/180);
b=0;
s0 = [a; b];
[t,s] = ode45(@pend_P,tspan,s0);
S1 = s(:,1);
S2 = s(:,2);
plot(t,S1*180/pi)


function sdot = pend_P(t,s)
a=10*(pi/180);
b=0;
% Expansion of s(1) or p :
p1 =0.000000001906*a^3 + (-0.0000007948)*a^2 + 0.00009188*a + (-0.003481);
p2 =0.00000915*a^2 + (-0.0009381)*a + 0.05331;
p3 =1*((-0.0001542)*a^2 + (-0.006078)*a + (-2.089));
p4 =a; 

s(1) =( (p1*t.^3) + (p2*t.^2) + (p3*t) + (p4) );

sdot = [s(2);  (2*s(2).^2/s(1))];     
end 

The curve I'm getting is a horixontal line though, as the value of S1 appears to be constant. Any idea where I've gone wrong?

Answer 1

[Alan Stevens]

Within your function pend you have

function sdot = pend_P(t,s)
a=10*(pi/180);
b=0;
...etc

This resets a and b to the initial conditions every time pend is called. You should have

function sdot = pend_P(t,s)
a=s(1);
b=s(2);
...etc

Also, if you start b (s(2)) at zero, then there is no way for anything to change given that

sdot = [s(2);  (2*s(2).^2/s(1))]; 

i.e. both terms in sdot will stay at zero.

[Nivedita Tanksali]

I tried what you suggested, but it made no difference.

Also, 'a' is the initial displacement, it is an initial condition. It is not equal to s(1). And b is not equal to s(2) either, as it is the inital velocity.

[Alan Stevens]

You pass a and b to the function as initial conditions, so, because b is zero, this means s(2) is intially zero, which means sdot returns 0, which means nothing changes! Because of the structure of sdot, the initial value of s(2) must be non-zero if you want anything to change.

If you intend that the values of p are calculated by the initial value of a only, and are unchanged after that, they are better calcuated outside of pend, and then passed in, otherwise they are recaculated on every call to pend, which is inefficient.