Add polynomial trendline to plot using script (m file) and code optimisation?

Hi,

I'd like to be able to add a trendline to my plot on figure(2) coding it in the m file. How would I go about doing that? It's currently only coming up as data points. My code is included below. Also, any potential optimizations in what I've currently written?

% Propulsion & Turbomachinery Lab
  %%Define Variables
length_moment_arm = 94; % mm
pressure = 1.02; % bar
density = (1.02*10^5)/(287*(19+273)); %kg/m^3
rotate_speed_n_rev = 100; % rev/s
rotate_speed_n_rad = rotate_speed_n_rev*2*pi; % rad/s
%%Thrust Calibration
thrust_temp = 19; % degrees
thrust_volts = [-0.1302, -0.3611, -0.2429, -0.8304, -0.482,...
      -0.7148, -0.9475, -0.1875, -1.0636, -0.5989]';
calibration_mass_g = [0, 100, 50, 300, 150, 250, 350, 25, ...
      400, 200];
calibration_mass_kg = calibration_mass_g(:)/1000;
calibration_force = calibration_mass_kg(:)*9.81;
thrust_volts = sort(thrust_volts,1,'descend');
calibration_force = sort(calibration_force);
sens_thrust = ((thrust_volts(10)-thrust_volts(1))/...
      (calibration_force(10)-calibration_force(1)));
c_thrust = thrust_volts(1);
%%Torque Calibration
torque_temp = 18; % degrees
torque_volts = [0.1180, 0.5823, 0.3506, 1.5138, 0.8185, 1.2833,...
      1.9809, 1.7476, 0.234, 1.0510];
torque_volts = sort(torque_volts);
sens_torque = ((torque_volts(10)-torque_volts(1))/...
      (calibration_force(10)-calibration_force(1)));
c_torque = torque_volts(1);
%%4 inch (10cm) pitch
prop1_d_cm = 25; % cm
prop1_d_m = prop1_d_cm/100; % m
prop1_thrust_offset_start = -0.0717; % Volts
prop1_thrust_offset_end = -0.1166; % Volts
prop1_thrust_offset_avg = (prop1_thrust_offset_start+prop1_thrust_offset_end)/2;
prop1_torque_offset_start = 0.0342; % Volts
prop1_torque_offset_end = 0.0353; % Volts
prop1_torque_offset_avg = (prop1_torque_offset_start+prop1_torque_offset_end)/2;
prop1_tunnel_velocity_min = 0:200:1200; % rev/min
prop1_tunnel_velocity_sec = prop1_tunnel_velocity_min/60; % rev/sec
prop1_thrust_volts = [-1.0438, -0.9779, -0.8025, -0.5777, -0.347, -0.1337, 0.1423]; % Volts
prop1_torque_volts = [0.331, 0.3344, 0.3103, 0.2617, 0.1975, 0.1246, 0.0176]; % Volts
prop1_delta_p = [0, 5, 16, 33, 60, 88, 137]; % Pa
prop1_tunnel_velocity_ms = (2.2*prop1_delta_p)/density).^0.5;
prop1_thrust = [prop1_thrust_volts-prop1_thrust_offset_avg]/sens_thrust; % Newtons
prop1_torque = [prop1_torque_volts-prop1_torque_offset_avg]/sens_torque; % Newtons
prop1_adv_ratio_J = prop1_tunnel_velocity_ms/(rotate_speed_n_rev*prop1_d_m);
prop1_CT = prop1_thrust/(density*(rotate_speed_n_rev^2)*(prop1_d_m^4));
prop1_CQ = prop1_torque/(density*(rotate_speed_n_rev^2)*(prop1_d_m^5));
prop1_power = prop1_torque*rotate_speed_n_rev;
prop1_prop_eff = (prop1_thrust.*prop1_tunnel_velocity_sec)./prop1_power;
%%Plotting section
figure(1)
subplot(1,2,1)
plot(calibration_force,thrust_volts)
grid on
title('Thrust Calibration')
xlabel('Force Applied (Newtons, N)')
ylabel('Thrust (Volts, V)')
subplot(1,2,2)
plot(calibration_force,torque_volts)
grid on 
title('Torque Calibration')
xlabel('Force Applied (Newtons, N)')
ylabel('Torque (Volts, V)')

figure(2)
for x = 1:6
  hold on
  plot(prop1_adv_ratio_J(1,x),prop1_prop_eff(1,x),'bo')
  title('Propulsive Effiency vs. Advance Ratio')
  xlabel('Advance Ratio, J')
  ylabel('Propulsive Effiency')
end

[dpb]

Not sure where you want; you mean a fitted curve in Fig 2? If so, use polyfit and polyval (probably with a quadratic it would appear).

For some simplification, no need for the loop at the end, simply write

plot(prop1_adv_ratio_J(1,1:6),prop1_prop_eff(1,1:6),'bo')
title('Propulsive Effiency vs. Advance Ratio')
xlabel('Advance Ratio, J')
ylabel('Propulsive Effiency')

[J]

Thanks!

So, having watched this video... https://www.youtube.com/watch?v=K2dsK3FFbik

I came up with this:

p1 = polyfit(prop1_adv_ratio_J(1,1:6),prop1_prop_eff(1,1:6),5)
t_interest = linspace(0,prop1_adv_ratio_J(6),100)
y_interest = polyval(p1,t_interest)
plot(t_interest,y_interest)
grid on

Works great, but how can I add now add the equation of the best fit line to a specific place on the plot? e.g. y = x^5 + x^4 + x^3 ...etc.

[dpb]

You used the fit from

p1=polyfit(J(1:6),eff(1:6),5);

(with variables/indices shortened to key differences for clarity...)

BAD IDEA!!! A fifth-order polynomial with only six datapoints is likely to have undesired inflection points in general. Turns out with your data it's not as bad as one might fear but try

x=linspace(0,prop1_adv_ratio_J(1,6),100);
p4=polyfit(prop1_adv_ratio_J(1,1:6),prop1_prop_eff(1,1:6),4);
y4=polyval(p4,x);
hold on
plot(x,y4,'r:')
p3=polyfit(prop1_adv_ratio_J(1,1:6),prop1_prop_eff(1,1:6),3);
y3=polyval(p4,x);
plot(x,y3,'g:')

The lesser orders have a little more prediction error but not the inflection of the 5th order between the first two points.

If you want an interpolating polynomial instead of a model, I suggest using an interpolating spline instead of one polynomial instead.

If

Answer 1

[Star Strider]

‘how can I add now add the equation of the best fit line to a specific place on the plot? e.g. y = x^5 + x^4 + x^3 ...etc.’

Use the text function or its friends.