ODE45 to solve vector ode

79 views (last 30 days)
Adam Binder
Adam Binder on 25 Jan 2020
Answered: ASHOK BANERJEE on 27 Oct 2020
I'm trying to solve the following differential equation
where I have and ,3x1 vectors, as intital conditions.
Here is my code:
clear all
clc
mu = 400000; %km^3/s^2
r0 = [9000;1400;800];%km
v0 = [-1.3;6.3;3.7];%km/s
tspan = [0;30000];
ic = [r0;v0];
f = @(t,y)[y(2);-(mu*y(1))/(norm(y(1))^3)];
[ts,ys] = ode45(f,tspan,ic);
I'm getting the error of:
@(T,Y)[Y(2);-(MU*Y(1))/(NORM(Y(1))^3)] returns a vector of length 2, but the length of initial conditions vector is 6. The
vector returned by @(T,Y)[Y(2);-(MU*Y(1))/(NORM(Y(1))^3)] and the initial conditions vector must have the same number of
elements.
Error in ode45 (line 115)
odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
Error in orbitplot (line 24)
[ts,ys] = ode45(f,tspan,ic);
How am I able to set up ode45 to be able to accept vectors as initial conditions?

Accepted Answer

James Tursa
James Tursa on 25 Jan 2020
Edited: James Tursa on 25 Jan 2020
You have a six element state vector. The y(1) and y(2) coming into your derivative function are not position and velocity vectors, they are the x and y position elements. It would probably be easier to write a short function for the derivative because of the calculations involved:
function dy = gravity(y,mu)
R = y(1:3);
V = y(4:6);
Rmag = norm(R);
RDOT = V;
VDOT = -mu * R / Rmag^3;
dy = [RDOT;VDOT];
end
Then call ode45 with a function handle that passes y and mu to the gravity function:
[ts,ys] = ode45(@(t,y)gravity(y,mu),tspan,ic);
ode45 isn't well suited for the orbit problem. You may have to use odeset( ) to create some tight tolerances to pass in to ode45 to get good results.
  2 Comments
Adam Binder
Adam Binder on 25 Jan 2020
Edited: Adam Binder on 26 Jan 2020
Thank you! So I need to have y(1:n) within the state space model for n initial condition elements?
James Tursa
James Tursa on 25 Jan 2020
Yes. The value of n depends on the number of equations and the order of the equations.

Sign in to comment.

More Answers (1)

ASHOK BANERJEE
ASHOK BANERJEE on 27 Oct 2020
[t,x]=ode45(@func,[0 5],[5 3]);
plot(t,x(:,1))
xlabel('Time[sec]')
ylabel('x(t)')
-----------------------------------------------------------------------------
function xdot = func(t,x)
%UNTITLED2 Summary of this function goes here
% Detailed explanation goes here
xdot(1)=x(2);
xdot(2)=20-7*x(2)-10*x(1);
xdot=xdot*;
end
Error in ode45 (line 115)
odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
-----------------------------------------------------------------------------

Products


Release

R2018b

Community Treasure Hunt

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

Start Hunting!