Symbolic Paramaters Not Supported in Nonpolynomial Equations

I am trying to use Matlab to solve a system of equations to determine the forces and moments acting on a satellite.

The variables F_d, q and q_dot are a vector of values, that represent force, pitch and pitch rate over time. Ideally, I would like this code to start with a value of theta to represent the starting pitch, and using the solutions of these equations model the change in pitch over time.

I am trying to solve the system of equations 1 - 4, at each time step t. However vpasolve is not working for me and I am struggling to think of what else to do? Any suggestions?

 for t = 1:100
    syms U_dot U W W_dot M_a theta
    % Force in X direction
    eqn1 = [-F_d(t) *cos(theta) - m_sat*g *sin(theta) == m_sat*(U_dot + q(t)*W)]; 
    %Force in Z direction
    eqn2 = [F_d(t)*sin(theta) + m_sat*g*cos(theta) == m_sat*(W_dot - q(t)*U)]; 
    %Pitching Moment
    eqn3 = [M_a == Iyy*q_dot(t)]; 
    eqn4 = [U == v_inf/cos(theta)];

    P = vpasolve([eqn1 eqn2 eqn3 eqn4],[U_dot W W_dot M_a theta]);
end 

Answer 1

[Torsten]

Use an ODE integrator (e.g. ODE45) to solve

dU/dt = -Fd(t)/m_sat * v_inf/U - g*sqrt(1-(v_inf/U)^2) - q(t)*W            (1')
dW/dt =  Fd(t)/m_sat * sqrt(1-(v_inf/U)^2) + g*v_inf/U + q(t)*U            (2')

with appropriate initial conditions for U and W.

These equations follow from the relation

cos(theta) = v_inf/U

It follows

sin(theta) = +/- sqrt(1-cos^2(theta)) = +/- sqrt(1-(v_inf/U)^2)

I used the positive root for the equations above:

sin(theta) = + sqrt(1-(v_inf/U)^2)

Equation 3 is superfluous.

You can use your arrays Ft(t) and q(t) to interpolate the values to the times requested by the integrator.

Answer 2

[Star Strider]

Solve it once, generate numeric functions, then evaluate the functions.

%  for t = 1:100
    syms U_dot U W W_dot M_a theta 
    syms F_d(t) m_sat g q(t) Iyy q_dot(t) v_inf 
    % Force in X direction
    eqn1 = (-F_d(t) *cos(theta) - m_sat*g *sin(theta) == m_sat*(U_dot + q(t)*W)); 
    %Force in Z direction
    eqn2 = (F_d(t)*sin(theta) + m_sat*g*cos(theta) == m_sat*(W_dot - q(t)*U)); 
    %Pitching Moment
    eqn3 = (M_a == Iyy*q_dot(t)); 
    eqn4 = (U == v_inf/cos(theta));

    P = solve([eqn1 eqn2 eqn3 eqn4],[U_dot W W_dot M_a theta]);
    U_dot_fcn = matlabFunction(P.U_dot)

Warning: Function 'F_d' not verified to be a valid MATLAB function.

U_dot_fcn = function_handle with value:
    @(U,g,m_sat,t,v_inf)[-((v_inf.*F_d(t))./U-g.*m_sat.*sqrt(-1.0./U.^2.*v_inf.^2+1.0))./m_sat;-((v_inf.*F_d(t))./U+g.*m_sat.*sqrt(-1.0./U.^2.*v_inf.^2+1.0))./m_sat]

    W_fcn = matlabFunction(P.W)

W_fcn = function_handle with value:
    @()[0.0;0.0]

    W_dot_fcn = matlabFunction(P.W_dot)

Warning: Function 'F_d' not verified to be a valid MATLAB function.

Warning: Function 'q' not verified to be a valid MATLAB function.

W_dot_fcn = function_handle with value:
    @(U,g,m_sat,t,v_inf)[(-U.*F_d(t).*sqrt(1.0./U.^2.*(U.^2-v_inf.^2))+U.^2.*m_sat.*q(t)+g.*m_sat.*v_inf)./(U.*m_sat);(U.*F_d(t).*sqrt(1.0./U.^2.*(U.^2-v_inf.^2))+U.^2.*m_sat.*q(t)+g.*m_sat.*v_inf)./(U.*m_sat)]

    M_a_fcn = matlabFunction(P.M_a)

Warning: Function 'q_dot' not verified to be a valid MATLAB function.

M_a_fcn = function_handle with value:
    @(Iyy,t)[Iyy.*q_dot(t);Iyy.*q_dot(t)]

    theta_fcn = matlabFunction(P.theta)

theta_fcn = function_handle with value:
    @(U,v_inf)[-acos(v_inf./U);acos(v_inf./U)]

% end 

.

[Rhiannon Elliott]

Hi, I had a go at this method but it didnt really work. I was not able to extract any values from it. I would still consider myself a beginner a MATLAB so I apologise for asking silly questions! I have included more of my code to hopefully provide better context.

close all
clear 
clc

%Aim of this code is to see if the thingy will converge and give me values
%with a range of density data changing over time.

%% Inputs
%Dimensions
x_sat = 1;                  %Satellite Width [m]
y_sat = 1;                  %Satellite Depth [m]
z_sat = 1;                  %Satellite Height[m]
m_sat = 200;                  %Satellite Mass [kg]

z_sail = 15;                 %Drag Sail Width [m]
y_sail = 15;                 %Drag Sail Depth [m]

m_sail = 2;                 %Drag Sail Mass [kg]

l_boom = 10;                 %Boom length [m]

Ixx = (1/12) * m_sat * (y_sat^2 + z_sat^2) + (1/12)*m_sail*(y_sail^2 + z_sail^2);
Iyy = (1/12) * m_sat * (x_sat^2 + z_sat^2) + ((1/12)*m_sail*z_sail^2) + (m_sail*l_boom^2);
Izz = (1/12) * m_sat * (x_sat^2 + y_sat^2) + ((1/12)*m_sail*y_sail^2) + (m_sail*l_boom^2);

%Constants
R_e = 6400;                                     %Radius of Earth [km]
w_e = 7.27e-5;                                  %Angular Rotation of Earth [rad/s]
mu = 3.98e14;                                   %Gravitational paramater of Earth [m^3/s^2]
Alt = [300];                %Altitudes of interest [km]
r = (Alt + R_e)*1000;                           %Radius [m]
A = z_sail*y_sail;                              %Area [m^2]
C_d = 1.1;                                      %Drag Coefficient for a flat plate [-]
theta_i = (10*pi)/180;                          %Initial angle of rotation of drag sail [rads]
g = 9.81;
%% Calculations
v_orbital = sqrt(mu./r);                         %Orbital Velocity [m/s]
v_inf = v_orbital - (w_e .*r);                   %V_infinity [m/s]

%% Density/Force calculations
%Solar Flux Values
%All values
%F10_7_daily, F10_7Avg, AP daily, ap_00_03_hr_prior, ap_03_06_hr_prior, ap_06_09_hr_prior, ap_09_12_hr_prior, ap_09_12_hr_prior, ap_33_59_hr_prior
SF = [78.5   79.8    2.5    2.0    2.0    3.0    2.0    2.1    7.0];
Period_100 = linspace(1,100,100);               %100s
[~,rho_300] = atmosnrlmsise00(300000, -55, 45 , 2021,1,(43200 + Period_100), SF(1), SF(2),SF(3:9));
rho_300 = rho_300(:,6);
F_d300 = 0.5*C_d*rho_300*A*abs(v_inf)*v_inf;
F_d = F_d300;

t = Period_100';
q = sqrt(F_d./(l_boom*m_sail));
dq = q(2:100) - q(1:99);
dt = t(2:100) - t(1:99);
q_dot = dq./dt;
q_dot = [0;q_dot];



for t = 1:100
    syms U_dot U W W_dot M_a theta 
    syms F_d(t) m_sat g q(t) Iyy q_dot(t) v_inf 
    % Force in X direction
    eqn1 = (-F_d(t) *cos(theta) - m_sat*g *sin(theta) == m_sat*(U_dot + q(t)*W)); 
    %Force in Z direction
    eqn2 = (F_d(t)*sin(theta) + m_sat*g*cos(theta) == m_sat*(W_dot - q(t)*U)); 
    %Pitching Moment
    eqn3 = (M_a == Iyy*q_dot(t)); 
    eqn4 = (U == v_inf/cos(theta));
    P = solve([eqn1 eqn2 eqn3 eqn4],[U_dot W W_dot M_a theta]);
   
    U_dot_fcn = matlabFunction(P.U_dot)
    W_fcn = matlabFunction(P.W)
    W_dot_fcn = matlabFunction(P.W_dot)
    M_a_fcn = matlabFunction(P.M_a)
    theta_fcn = matlabFunction(P.theta)
end

[Star Strider]

One problem is that ‘U’ is not assigned any value, so all calculations using it will fail. Provide a scalar or appropriately-sized vector for it, and the functions should produce results.

Also, once ‘U’ has a value, replace the vpa calls with double to produce values that can be used in other parts of the code, or saved to an appropriate file.

%Aim of this code is to see if the thingy will converge and give me values
%with a range of density data changing over time.

% for t = 1:100                             % Create The Functions First
sympref('AbbreviateOutput',false);
    syms U_dot U W W_dot M_a theta 
    syms F_d m_sat g q Iyy q_dot v_inf 
    % Force in X direction
    eqn1 = (-F_d *cos(theta) - m_sat*g *sin(theta) == m_sat*(U_dot + q*W)); 
    %Force in Z direction
    eqn2 = (F_d*sin(theta) + m_sat*g*cos(theta) == m_sat*(W_dot - q*U)); 
    %Pitching Moment
    eqn3 = (M_a == Iyy*q_dot); 
    eqn4 = (U == v_inf/cos(theta));
    P = solve([eqn1 eqn2 eqn3 eqn4],[U_dot W W_dot M_a theta]);
    
    U_dot_fcn = matlabFunction(P.U_dot)

U_dot_fcn = function_handle with value:
    @(F_d,U,g,m_sat,v_inf)[-((F_d.*v_inf)./U+g.*m_sat.*sqrt(-1.0./U.^2.*v_inf.^2+1.0))./m_sat;-((F_d.*v_inf)./U-g.*m_sat.*sqrt(-1.0./U.^2.*v_inf.^2+1.0))./m_sat]

%     W_fcn = matlabFunction(P.W)
    W_dot_fcn = matlabFunction(P.W_dot)

W_dot_fcn = function_handle with value:
    @(F_d,U,g,m_sat,q,v_inf)[(g.*m_sat.*v_inf+F_d.*U.*sqrt(1.0./U.^2.*(U.^2-v_inf.^2))+U.^2.*m_sat.*q)./(U.*m_sat);(g.*m_sat.*v_inf-F_d.*U.*sqrt(1.0./U.^2.*(U.^2-v_inf.^2))+U.^2.*m_sat.*q)./(U.*m_sat)]

    M_a_fcn = matlabFunction(P.M_a)

M_a_fcn = function_handle with value:
    @(Iyy,q_dot)[Iyy.*q_dot;Iyy.*q_dot]

    theta_fcn = matlabFunction(P.theta)

theta_fcn = function_handle with value:
    @(U,v_inf)[acos(v_inf./U);-acos(v_inf./U)]

% end

%% Inputs
%Dimensions
x_sat = 1;                  %Satellite Width [m]
y_sat = 1;                  %Satellite Depth [m]
z_sat = 1;                  %Satellite Height[m]
m_sat = 200;                  %Satellite Mass [kg]

z_sail = 15;                 %Drag Sail Width [m]
y_sail = 15;                 %Drag Sail Depth [m]

m_sail = 2;                 %Drag Sail Mass [kg]

l_boom = 10;                 %Boom length [m]

Ixx = (1/12) * m_sat * (y_sat^2 + z_sat^2) + (1/12)*m_sail*(y_sail^2 + z_sail^2);
Iyy = (1/12) * m_sat * (x_sat^2 + z_sat^2) + ((1/12)*m_sail*z_sail^2) + (m_sail*l_boom^2);
Izz = (1/12) * m_sat * (x_sat^2 + y_sat^2) + ((1/12)*m_sail*y_sail^2) + (m_sail*l_boom^2);

%Constants
R_e = 6400;                                     %Radius of Earth [km]
w_e = 7.27e-5;                                  %Angular Rotation of Earth [rad/s]
mu = 3.98e14;                                   %Gravitational paramater of Earth [m^3/s^2]
Alt = [300];                %Altitudes of interest [km]
r = (Alt + R_e)*1000;                           %Radius [m]
A = z_sail*y_sail;                              %Area [m^2]
C_d = 1.1;                                      %Drag Coefficient for a flat plate [-]
theta_i = (10*pi)/180;                          %Initial angle of rotation of drag sail [rads]
g = 9.81;
%% Calculations
v_orbital = sqrt(mu./r);                         %Orbital Velocity [m/s]
v_inf = v_orbital - (w_e .*r);                   %V_infinity [m/s]

%% Density/Force calculations
%Solar Flux Values
%All values
%F10_7_daily, F10_7Avg, AP daily, ap_00_03_hr_prior, ap_03_06_hr_prior, ap_06_09_hr_prior, ap_09_12_hr_prior, ap_09_12_hr_prior, ap_33_59_hr_prior
SF = [78.5   79.8    2.5    2.0    2.0    3.0    2.0    2.1    7.0];
Period_100 = linspace(1,100,100);               %100s
[~,rho_300] = atmosnrlmsise00(300000, -55, 45 , 2021,1,(43200 + Period_100), SF(1), SF(2),SF(3:9));
rho_300 = rho_300(:,6);
F_d300 = 0.5*C_d*rho_300*A*abs(v_inf)*v_inf;
F_d = F_d300;

t = Period_100';
q = sqrt(F_d./(l_boom*m_sail));
dq = q(2:100) - q(1:99);
dt = t(2:100) - t(1:99);
q_dot = dq./dt;
q_dot = [0;q_dot];

U_dot_v = vpa(U_dot_fcn(F_d,U,g,m_sat,v_inf))

U_dot_v = 

W_dotv = vpa(W_dot_fcn(F_d,U,g,m_sat,q,v_inf))

W_dotv = 

M_av = vpa(M_a_fcn(Iyy,q_dot))

M_av = 

theta_v = vpa(theta_fcn(U,v_inf))

theta_v = 

.

[Torsten]

@Star Strider

The OP defines two ODEs for U and W with time-dependent inputs F_d and q.

See my answer above.