Reputation: 105
So I need to solve x''(t) = -x(t)^p with initial conditions x(0)= 0 and v(0) = x'(0) = v_o = 1. The value of the parameter p is 1.
This is what I have:
function [t, velocity, x] = ode_oscilation(p)
y=[0;0;0];
% transform system to the canonical form
function y = oscilation_equation(x,p)
y=zeros(2,1);
y(1)=y(2);
y(2)=-(x)^p;
% to make matlab happy we need to return a column vector
% so we transpose (note the dot in .')
y=y.';
end
tspan=[0, 30]; % time interval of interest
[t,velocity,x] = ode45(@oscilation_equation, tspan, 1);
t = y(:,1);
xposition=y(:,3);
velocity=y(:,2);
end
and this is the error message I receive:
ode_oscillation(1) Error using odearguments (line 91) ODE_OSCILLATION/OSCILATION_EQUATION must return a column vector.
Error in ode45 (line 114) [neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...
Error in ode_oscillation (line 17) [t,velocity,x] = ode45(@oscilation_equation, tspan,1);
Upvotes: 0
Views: 4403
Reputation: 38032
There's a few things going wrong here. First, from help ode45
:
ode45 Solve non-stiff differential equations, medium order method.
[TOUT,YOUT] = ode45(ODEFUN,TSPAN,Y0) with TSPAN = [T0 TFINAL] integrates the system of differential equations y' = f(t,y) from time T0 to TFINAL with initial conditions Y0.
Note that ode45
expects a function f(t,y)
, where size(t) == [1 1]
for time and size(y) == [1 N]
or [N 1]
for solution values. Your oscilation_equation
has the order of input arguments inverted, and you input a constant parameter p
instead of time t
.
Also, the initial conditions Y0
should have the same size as y
; so size(y0) == [N 1]
or [1 N]
. You just have 1
, which is clearly causing errors.
Also, your output arguments t
, xposition
and velocity
will be completely ignored and erroneous, since y
is not set as output argument from ode45
, and most of all, their names do not correspond to ode_oscilation
's output arguments. Also, their order of extracting from columns of y
is incorrect.
So, in summary, change everything to this:
function [t, v, x] = ode_oscilation(p)
% initial values
y0 = [0 1];
% time interval of interest
tspan =[0 30];
% solve system
[t,y] = ode45(@(t,y) [y(2); -y(1)^p], tspan, y0);
% and return values of interest
x = y(:,1);
v = y(:,2);
end
Upvotes: 3