Error trying to solve 2 Second ODE

1 view (last 30 days)
Scott Angus
Scott Angus on 15 Nov 2019
Commented: Scott Angus on 15 Nov 2019
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!
  1 Comment
darova
darova on 15 Nov 2019
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;

Sign in to comment.

Answers (1)

James Tursa
James Tursa on 15 Nov 2019
Edited: James Tursa on 15 Nov 2019
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) ];
  2 Comments
Scott Angus
Scott Angus on 15 Nov 2019
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
Scott Angus on 15 Nov 2019
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.

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!