0

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

Oliver Amundsen
  • 1,491
  • 2
  • 21
  • 40
  • 1
    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 Feb 05 '13 at 23:56
  • 1
    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 Feb 06 '13 at 00:03
  • By the way: It is not best practice to use global variables here. For transfer of Parameters with ode45 see: http://stackoverflow.com/questions/29215121/matlab-for-loop-within-ode-solver/29219310#29219310 – BerndGit Apr 03 '15 at 12:11
  • thanks pal , that's certainly helpful ! – Oliver Amundsen Apr 03 '15 at 15:17

2 Answers2

2

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.

s.bandara
  • 5,636
  • 1
  • 21
  • 36
  • You were right they were empty. Now it works and I get the value of both variables, but I get this error message: ??? Error using ==> odearguments at 113 YPRIME must return a column vector.... – Oliver Amundsen Feb 06 '13 at 00:21
0

the syntax for ODE solvers is [s y]=ode45(@yprime, [1 10], [2 2])

and you dont need to do elementwise operation in your case i.e. instead of .* just use *y(:,1) plotted against s with initial conditions [2 2] and parameter values: g=1;a=2;v0=1;d=1.5;

y(:,2) vs s

WYSIWYG
  • 494
  • 6
  • 23
  • i tried and it just worked.. just didnt know suitable values for your global variables... it works fine.. – WYSIWYG Feb 08 '13 at 04:23
  • well in your equation there is a division by s; so you cannot start s from 0. just missed it.. and perhaps thats why you are getting NaN or something – WYSIWYG Feb 08 '13 at 05:04
  • Thanks bharat_iyengar. The values of the global variables are vectors of 1x100, so that's why I used the .*, otherwise it will say: "??? Error using ==> mpower Inputs must be a scalar and a square matrix" – Oliver Amundsen Feb 11 '13 at 17:06