Error in ode45

7 vues (au cours des 30 derniers jours)
christian fares
christian fares le 18 Mai 2020
I added to a code another planet and when i added the Mercury planet on ode45 i got an error that is:
Error in ode45 (line 115)
odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
Error in SimpleGravity (line 39)
[attime,finalposition]=ode45(@Propagator,[0:stepsize:finaltime],[EarthInitial]);
%% Inputs
stepsize=1; %step size in days, should not exceed 1
finaltime=10; %days from t0
%% System Conditions
global sunm G
sunm=1.9891e30; %kg
G=1.4879*10^(-34); %AU^3/(kg*Day^2)
%% Earth Initial Conditions
x=-7.829690890802676E-01; %AU
y=6.009223815219861E-01; %AU
z=-1.377986550047871E-05; %AU
Vx=-1.075646922487205E-02;%AU/day
Vy=-1.370584048238274E-02;%AU/day
Vz=3.974096024543801E-08; %AU/day
EarthInitial=[x,y,z,Vx,Vy,Vz]; %summary of initial conditions of the earth
%% New planet
Mx=-2.829690890802676E-01; %AU
My=8.009223815219861E-01; %AU
Mz=-4.377986550047871E-05; %AU
MVx=-0.075646922487205E-02;%AU/day
MVy=-2.370584048238274E-02;%AU/day
MVz=5.974096024543801E-08; %AU/day
MercuryInitial=[Mx,My,Mz,MVx,MVy,MVz];
%% Call the Function
%outpu=call the function, time to run, initial conditions
[attime,finalposition]=ode45(@Propagator,[0:stepsize:finaltime],[EarthInitial,MercuryInitial]); %% this is where the problem is when I added Mercuryinitial
%% Plot the Result
hold ON
plot3(finalposition(:,4),finalposition(:,5),finalposition(:,6),'g'); %plot of position of earth in green
plot3(finalposition(:,10),finalposition(:,11),finalposition(:,12),'g');
xlabel('Distance (AU)');ylabel('Distance (AU)');zlabel('Distance (AU)'); %label the axis
legend('Earth'); %tell me what the green line is
axis equal %fix the axis so it has the same scale in x, y and z
%note: sun is at 0,0,0
  2 commentaires
Steven Lord
Steven Lord le 18 Mai 2020
That doesn't look like the whole error message. Please show us all the text displayed in red and/or orange when you try to run your code.
christian fares
christian fares le 18 Mai 2020
Modifié(e) : christian fares le 18 Mai 2020
Oh im sorry your right! this is everything and thank you for your help
This is all the error:
Error using *
Incorrect dimensions for matrix multiplication. Check that the number of columns in the first matrix matches the number of rows in the second matrix. To perform
elementwise multiplication, use '.*'.
Error in Propagator (line 37)
TotAccelEarth=A1*input+A2; %sum up the linear + nonlinear accelerations
Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.
Error in ode45 (line 115)
odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
Error in SimpleGravity (line 39)
[attime,finalposition]=ode45(@Propagator,[0:stepsize:finaltime],[EarthInitial,MercuryInitial ]);
This is the function propagator:
function [out]=Propagator(t,input);
global sunm G
%state space representation idot=A1*input+A2(input)
%where idot is the derivative of the input
%this represents the linear part in A1 added to the non-linear part in A2
% |Vx| |0 0 0 1 0 0| |x | | 0 |
% |Vy| |0 0 0 0 1 0| |y | | 0 |
% |Vz| = |0 0 0 0 0 1| * |z | + | 0 |
% |Ax| |0 0 0 0 0 0| |Vx| | -G*sunm/R^2*|Rx| |
% |Ay| |0 0 0 0 0 0| |Vy| | -G*sunm/R^2*|Ry| |
% |Az| |0 0 0 0 0 0| |Vz| | -G*sunm/R^2*|Rz| |
% idot = A1 * input + A2(input)
% where |Rx| is the unit vector of R in the x direction
% where R=sqrt(x^2+y^2+z^2);
dx=zeros(6,1); %initialize the xdot matrix
A1=zeros(6,6);
A2=zeros(6,1);
A1(1,4)=1;
A1(2,5)=1;
A1(3,6)=1;
%% Earth
x=input(1); %x position
y=input(2); %y position
z=input(3); %z position
R = norm([x,y,z]); %find the radius from the earth to the sun
unitvector=[x,y,z]/R;
% do these calculations now so we only have to calculate it once
A2(4)=-G*sunm/R^2*unitvector(1); %acceleration in the x
A2(5)=-G*sunm/R^2*unitvector(2); %acceleration in the y
A2(6)=-G*sunm/R^2*unitvector(3); %acceleration in the z
TotAccelEarth=A1*input+A2; %sum up the linear + nonlinear accelerations
out=TotAccelEarth; %output the accelerations
%% Mercury
Mx=input(7); %x position
My=input(8); %y position
Mz=input(9); %z position
R = norm([Mx,My,Mz]); %find the radius from the earth to the sun
unitvector=[Mx,My,Mz]/R;
A2(4)=-G*sunm/R^2*unitvector(1); %acceleration in the x
A2(5)=-G*sunm/R^2*unitvector(2); %acceleration in the y
A2(6)=-G*sunm/R^2*unitvector(3); %acceleration in the z
TotAccelMercury=A1*input+A2; %sum up the linear + nonlinear accelerations
out=TotAccelMercury; %output the accelerations
end

Connectez-vous pour commenter.

Réponse acceptée

Steven Lord
Steven Lord le 18 Mai 2020
Your initial condition vector has twelve elements, so the input to your ODE function Propagator will be a 12-by-1 vector. [As an aside, I would use a different name for the input argument than "input". For one thing it says effectively nothing about what that vector represents. Something like planetPositions would be more descriptive. For another, input already has a meaning in MATLAB. While you wouldn't want to call input in your ODE function because it would prompt the user at each timestep, it's a general good habit to get into not to use MATLAB function names for your variables.]
Your A1 variable is a 6-by-6 matrix. You can't multiply a 6-by-6 matrix and a 12-by-1 vector.
Perhaps you want to reshape the input to be a 6-by-2 matrix, perform the multiplication to receive a 6-by-2 matrix, and reshape the result back to 12-by-1 so it's the correct size for your ODE function to return?
  2 commentaires
christian fares
christian fares le 18 Mai 2020
Okay , this is what i changed: ( i changed the name to planetpositions and reshaped it to be a 6-by-2 matrix)
but i didnt know how to reshape the result back to 12-by-1.
function [out]=Propagator(t,PlanetPositions);
PlanetPositions= zeros(6,2);
Errors:
Error using odearguments (line 93)
PROPAGATOR must return a column vector.
Error in ode45 (line 115)
odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
Error in SimpleGravity (line 39)
[attime,finalposition]=ode45(@Propagator,[0:stepsize:finaltime],[EarthInitial,MercuryInitial ]);
Walter Roberson
Walter Roberson le 18 Mai 2020

Connectez-vous pour commenter.

Plus de réponses (1)

Meysam Mahooti
Meysam Mahooti le 28 Juil 2021
Modifié(e) : Walter Roberson le 29 Juil 2021

Catégories

En savoir plus sur Programming dans Help Center et File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by