I am unable to find out ans to this code please tell me how to remove these errors. errors are like: Error using sym>checkindex (line 1545) Index must be a positive integer or logical. Error in sym>privformatscalar (line 1492) checkindex(x); Er

clear all
n = 172
H = 20000
CF = 0.5
V_wind =10
%%Constants taken
u = [7.447323    2.071808    9.010095    7.981491 194];
a = u(1); % 7.447323;
b = u(2); % 2.071808;
c = u(3); % 9.010095;
d = u(4); % 7.981491;
l = u(5); % 194;
%%Profile
x = (0:0.1:u(5));
f = 1/8.*(sqrt(u(1).*(u(5)-x).*(u(2).*x-u(5)*sqrt(u(3))+sqrt(u(3)*u(5).^2-u(4)*u(5).*x))));
D = 2*max(f);
R = D/2;
clear x

%%Area of Surface
syms x 
y = 1/8*(sqrt(u(1)*(u(5)-x).*(u(2).*x-u(5)*sqrt(u(3))+sqrt(u(3)*u(5).^2-u(4)*u(5).*x))));   % Given curve for shape
y1 = diff(y);
y2 = sqrt(1+y1.^2);
y3 = y*y2;
A = 2*pi*int(y3,0,u(5));
vpa(A,5)                                                        % Area of Envelope
A = ans

m = 120*10.^-3*A*0.12*10.^-3

clear x y
%%constants
alpha = 0.15;
albedo_ground = 0.35;
albedo_cloud = 0.5;
P_sea = 101325;
e = 0.0167;
%%solar model
syms t
delta = 23.45.*sin((pi/180).*(360/365).*(284+n));
% [T_oa, P_oa, rho_oa]  = atmos(H);
P_oa = 5446;
rho_oa = 0.088;
T_oa = 216.65;
w = 15.*(12-t);
phi = 18.975;

theta = asin(sin((pi/180).*delta).*sin((pi/180).*phi)+cos((pi/180).*delta).*cos((pi/180).*phi).*cos((pi/180).*w));
% psi = asin((sin(pi/180).*w).*cos((pi/180).*delta)/cos(theta));
MA = 2.*pi.*n/365;
M = (P_oa/P_sea).*(sqrt(1229+(614.*sin(theta)).^2)-614.*sin(theta));  %call fun
Tau_atm = 0.5.*(exp(-0.65.*M)+exp(-0.95.*M));
zita = MA+2.*e.*sin(M)+1.25.*e.^2.*sin(2.*M); 
ID = (1367.*Tau_atm.^M).*((1+e.*cos(zita))./(1-e.^2)).^2;
I_d = (1/2).*ID.*sin(theta).*(1.-Tau_atm.^M)./(1.-1.4.*log(Tau_atm));
Q_D = alpha.*ID.*sin(theta)                                  % direct radiation

% diffuse soalr radiation model

Q_d = alpha.*I_d

%Reflected solar radiation model
albedo = albedo_ground.*(1-CF).^2+albedo_cloud.*CF;
Ir = albedo.*(ID.*sin(theta)+I_d);
omega = pi.*abs(12-t)/12; 
Q_r = alpha.*Ir.*cos(omega)

%Infrared Radiation from earth
T_g = 280;                                                 %temperature of ground in k
T_cl = 262.15;                                             % assuming cloud at altitude 4km in k
R = 6371000;                                               %redius of earth
epsilon_e = 0.95;      
alpha_ir = 0.85;
sigma = 5.67.*10.^-8;
tau_atm_IR = 1.716-0.5.*(exp(-0.65.*P_oa/P_sea) + exp(-0.95.*P_oa/P_sea));
T_e = T_g.*(1-CF)+T_cl.*CF;
eta = asin(R/(R+H));
VF = (1-cos(eta))/2;
Q_earth = sigma.*epsilon_e.*alpha_ir.*VF.*tau_atm_IR.*T_e.^4

%Infrared Radiation from sky
Pw = 0.0042;  
epsilon_sky = 1;
T_sky = T_oa.*(0.51+0.208.*sqrt(Pw)).^0.25;
Q_sky = sigma.*epsilon_sky.*alpha_ir.*(1-VF).*T_sky.^4

%Infrared Radiation from airship
syms T_f
epsilon_f = 0.85;
Q_a = 2.*sigma.*epsilon_f.*T_f.^4

r = 0.05;
Q_aIR = alpha_ir.*sigma.*epsilon_f.*T_f.^4/(1-r)

% Convection heat transfer model
% external convection
g = 9.81;
K_air = 0.0241.*(T_oa/273.15).^0.9;
beta_air = 1/T_oa;
Pr_air = 0.804-3.25.*10.^(-4).*T_oa;
nu_air= 1.46.*10.^-6.*T_oa.^1.5/(rho_oa.*(T_oa+110.4));
Ra_air = g.*beta_air.*(T_f-T_oa).*D.^3/nu_air.^2.*Pr_air;
Re_air = V_wind.*l/nu_air;
h_free = (K_air/D).*(0.6+0.387.*((Ra_air)/(1+(0.559/Pr_air).^(9/16)).^(16/9)).^(1/6)).^2;
h_forced = (K_air/l).*(Re_air.*Pr_air.*(0.2275/(log(Re_air)).^2.584-850/Re_air));
h_oc = (h_free.^3 + h_forced.^3).^(1/3);
Q_oc = h_oc.*(T_f-T_oa)

% %internal convection
% syms T_he
% rho_he = 0.1786
% k_he = 0.144.*(T_he/273.15).^0.7
% Pr_he = 0.729-1.6.*10.^(-4).*T_he
% beta_he = 1/T_he
% nu_he = 4.32.*10.^-6.*(T_he.^0.674)/rho_he
% Ra_he = g.*beta_he.*(T_f-T_he).*D.^3/nu_he.^2.*Pr_he
% if Ra_he<1.5.*10.^8
%     Nu_ic = 2.5.*(2+0.6.*Ra_he.^(1/4))
% else Nu_ic = 0.325.*Ra_he.^(1/3)
% end
%     
% h_ic = Nu_ic.*k_he/D 
% Q_ic=h_ic.*(T_f-T_he)

%%solution
C_p = 1.42;
h = 1;
x = 0:h:23
y = zeros(1,length(x))
y(1) = 216.65
F_xy = (Q_D+Q_d+Q_r+Q_earth+Q_sky-Q_aIR-Q_a-Q_oc)/(m.*C_p)

for i=1:(length(t))                              % calculation loop
    k_1 = F_xy(x(i),y(i))
    k_2 = F_xy(x(i)+0.5.*h,y(i)+0.5.*h.*k_1)
    k_3 = F_xy((x(i)+0.5.*h),(y(i)+0.5.*h.*k_2))
    k_4 = F_xy((x(i)+h),(y(i)+k_3.*h))

      y(i+1) = y(i) + (1/6).*(k_1+2.*k_2+2.*k_3+k_4).*h  % main equation
  end

Answer 1

[Walter Roberson]

You did not indicate which line the problem occurred on, or if you did then it got cut off because you put your entire question in the title.

Without running the program, I can predict that the difficult you have is that you have defined a symbolic formula, say f, and you are attempting to evaluate the symbolic formula with a specific numeric value, say x = 3.2343, by using the syntax f(x) . Symbolic formula are not functions, so MATLAB took that to be a request to index the variable at the given value.

If you want to evaluate a formula f written in terms of symbol x at a given numeric value, and get out a symbolic answer, then you should use

subs(f, 'x', the_numeric_value)

If you are expecting a symbolic number as a result and you want to get the double-precision equivalent, then

double( subs(f, 'x', the_numeric_value) )

If you want to evaluate a symbolic formula at a specific point and get out a double precision number and you are doing that repeatedly, then use matlabFunction() to convert the symbolic form to an anonymous function, and then evaluate that anonymous function at the numeric locations. For example,

Fn = matlabFunction(f);
Fn([7.2, 8.3])