how to solve a system of Ordinary Differential Equations (ODE's) in Matlab

Question

I have to solve a system of ordinary differential equations of the form:

dx/ds = 1/x * [y* (g + s/y) - a*x*f(x^2,y^2)]
dy/ds = 1/x * [-y * (b + y) * f()] - y/s - c

where x, and y are the variables I need to find out, and s is the independent variable; the rest are constants. I've tried to solve this with ode45 with no success so far:

y = ode45(@yprime, s, [1 1]);

function dyds = yprime(s,y)
    global g a v0 d
    dyds_1 = 1./y(1) .*(y(2) .* (g + s ./ y(2)) - a .* y(1) .* sqrt(y(1).^2 + (v0 + y(2)).^2));
    dyds_2 = - (y(2) .* (v0 + y(2)) .* sqrt(y(1).^2 + (v0 + y(2)).^2))./y(1) - y(2)./s - d;
   dyds = [dyds_1; dyds_2];
return

where @yprime has the system of equations. I get the following error message:

YPRIME returns a vector of length 0, but the length of initial conditions vector is 2. The vector returned by YPRIME and the initial conditions vector must have the same number of elements.

Any ideas? thanks

My bet would be that at least one of g, a, v0, or d remain uninitialized, thus []. Using these "coefficients" will yield an empty vector for dyds. You could test this with assert(~isempty(v0), 'v0 not initialized') in yprime. — s.bandara
Your equations are singular. You are dividing by x, but your initial condition is x = 0. I don't know if this is the source of your error, but it will be a problem. — Warren Weckesser
By the way: It is not best practice to use global variables here. For transfer of Parameters with ode45 see: stackoverflow.com/questions/29215121/… — BerndGit

s.bandara s.bandara · Accepted Answer · 2013-02-05T23:46:43

Certainly, you should have a look at your function yprime. Using some simple model that shares the number of differential state variables with your problem, have a look at this example.

function dyds = yprime(s, y)
    dyds = zeros(2, 1);
    dyds(1) = y(1) + y(2);
    dyds(2) = 0.5 * y(1);
end

yprime must return a column vector that holds the values of the two right hand sides. The input argument s can be ignored because your model is time-independent. The example you show is somewhat difficult in that it is not of the form dy/dt = f(t, y). You will have to rearrange your equations as a first step. It will help to rename x into y(1) and y into y(2).

Also, are you sure that your global g a v0 d are not empty? If any one of those variables remains uninitialized, you will be multiplying state variables with an empty matrix, eventually resulting in an empty vector dyds being returned. This can be tested with

assert(~isempty(v0), 'v0 not initialized');

in yprime, or you could employ a debugging breakpoint.

how to solve a system of Ordinary Differential Equations (ODE's) in Matlab

2 Answers