Error trying to solve 2 Second ODE

Hi,

I have a function that solves a differential equation using the RK method:

function [tvec, xvec] = rksolve(f, t0, tf, x0, dt)
%Solve f with initial conditions x0 between t0 -> tf
%with increments dt

%Set initial conditions
t = t0;
x = x0;

%Answer output
tvec = t;
xvec = x;


while t < tf
    %RK method of solving
    k1 = f(x, t);
    k2 = f(x+0.5*k1*dt, t+0.5*dt);
    k3 = f(x+0.5*k2*dt, t+0.5*dt);
    k4 = f(x+k3*dt, t+dt);
    x = x + ((k1 + 2*k2 + 2*k3 +k4)*dt)/6;
    t = t+dt;

    %Store outputs 
    tvec = [tvec t];
    xvec = [xvec x];
end

But when I try to run it with the equations defined in func.m, it comes up with a horzcat error in the line:

xvec = [xvec x];

func.m:

function a = coords(x, t);
%Initial Conditions
xpos = x(1);
vx = x(2);
ypos = x(3);
vy = x(4);

a = [(-xpos - 2*xpos*ypos); (-ypos - (xpos^2) +ypos^2)];

%dxdt = vx
%dvxdt = -x -2*x*y;

%dydt = vy;
%dvydt = -y -(x^2) + y^2;

Any help would be greatly appreciated!

[darova]

What about preallocation?

n = round((tf-t0)/dt+1);
tvec1 = zeros(1,n);
i = 1;
while t < tvec
    %% ...
    tvec1(i) = t;
    i = i + 1;
end
tvec = tvec1;

Answer 1

[James Tursa]

Just looking at func.m, it appears you pass in a 4-element x vector, but you only return a 2-element a vector. You need to return the derivative for all the states. So something like this instead based on your comments:

a = [ vx; (-xpos - 2*xpos*ypos); vy; (-ypos - (xpos^2) +ypos^2) ];

[Scott Angus]

Thanks for your response. I tried changing the output of func.m to

a = [ vx; (-xpos - 2*xpos*ypos); vy; (-ypos - (xpos^2) +ypos^2) ];

but still come up with an error:

>> [ts, xmat] = rksolve(@func, 0, 200, [0, 0.3, 0, 0], 0.2)
Error using horzcat
Dimensions of arrays being concatenated are not consistent.

Error in rksolve (line 25)
    xvec = [xvec x];

I'm not sure why because the x and xvec dimensions should be the same?

[Scott Angus]

I found the issue. I was passing the function rksolve with a horizontal vector for x0 instead of a vertical vector. Adding semicolons between the components of the vector when I passed it solved the problem.