ODE45 to solve vector ode

102 次查看（过去 30 天）
Adam Binder
Adam Binder 2020-1-25
回答： ASHOK BANERJEE 2020-10-27
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?

请先登录，再回答此问题。

采纳的回答

James Tursa
James Tursa 2020-1-25
编辑：James Tursa 2020-1-25
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 个评论
Adam Binder
Adam Binder 2020-1-25
编辑：Adam Binder 2020-1-26
Thank you! So I need to have y(1:n) within the state space model for n initial condition elements?
James Tursa
James Tursa 2020-1-25
Yes. The value of n depends on the number of equations and the order of the equations.

请先登录，再进行评论。

更多回答（1 个）

ASHOK BANERJEE
ASHOK BANERJEE 2020-10-27
[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);
-----------------------------------------------------------------------------

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

产品


版本

R2018b

Community Treasure Hunt

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

Start Hunting!

Translated by