ODE45 to solve vector ode
18 views (last 30 days)
Show older comments
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?
0 Comments
Accepted Answer
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
James Tursa
on 25 Jan 2020
Yes. The value of n depends on the number of equations and the order of the equations.
More Answers (1)
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);
-----------------------------------------------------------------------------
0 Comments
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!