Curve Fitting with Polynomials using Python

Topic: Method of Curve Fitting with Polynomial Functions;

Subject: Numerical Methods;

Tool: Python;

By: Gani Comia, Dec. 2024;

Curve fitting with polynomials finds a polynomial equation of a given degree that approximates the set of data points. Curve fitting is used in the following applications:

• Data Smoothing
• Predictive Modelling
• Scientific Analysis

The general form of a polynomial function of degree n is:

$$P(x)\;=\;{a_n\,x^n}+{a_{n-1}\,x^{n-1}}+\ldots+{a_1\,x}+{a_0}$$

This article shows a Python script as toolbox to implement the method of curve fitting with polynomials. The script requires numerical data for both dependent and independent variables. For a given sample of data points in Fig.1., an approximate polynomial curve of 3rd degree shown in Fig.2. can be calculated and plotted using Python script.

$$xdata = [ 0 \quad 1 \quad 2 \quad 3 \quad 4 \quad 5 \quad 6 \quad 7 ]$$

$$ydata = [ 5 \quad 10 \quad 60 \quad 50 \quad 130 \quad 220 \quad 300 \quad 280 ]$$

Fig.1. Sample Data Points for Approximating Polynomial Curve.

Below is the approximate polynomial curve that closely fit the given data points.

Fig.2. Polynomial Curve of 3rd Degree from Data Points.

Python Script

# -*- coding: utf-8 -*-
"""
Created on Mon Nov  7 21:08:48 2024
@author: Gani Comia
"""
'''
Curve Fitting with Polynomials for a given xdata and ydata
'''
import numpy as np
import matplotlib.pyplot as plt

# given: xdata and ydata
xdata = np.array([0. , 1. , 2. , 3. , 4. , 5. , 6. , 7.])
ydata = np.array([5. , 10. , 60. , 50. , 130. , 220. , 300. , 280.])

# solution:
# solving for coefficients of polynomial fit for cubic (degree=3)
# z is an array of coefficients , highest first , i.e.
# z = np.array ([c_3, c_2 , c_1 , c_0])
z = np.polyfit(xdata, ydata, 3)
print("\n",z,"\n") 

# ‘poly1d‘ returns the polynomial function from coefficients    
p = np.poly1d(z)
print(p)

# polynomial's curve data points
xs = [0.1 * i for i in range (70)]  # from 0 to 7, step=0.1
ys = [p(x) for x in xs]   # evaluate p(x) for all x in list xs

# creating plot of data points and polynomial curve for visualization
plt.plot(xdata, ydata, 'o', label='data points') # data points
plt.plot(xs, ys, label='fitted curve')    # polynomial fit
plt.title('Curve Fitting with Polynomials')
plt.ylabel('y-value')
plt.xlabel('x-value')
plt.legend()
plt.grid()
plt.show()

On the above script, Python Numpy’s utilities, np.polyfit(), returns the polynomial coefficients of the given degree that approximately fits on the given data points.  Numpy’s np.poly1d() creates a polynomial functions for plotting and calculation.

Feel free to comment for inquiry, clarification, or suggestion for improvement. Drop your email to request the softcopy of the file.

Disclaimer: The formulas and calculations presented are for technical reference only. Users must verify the accuracy and ensure compliance with applicable engineering standards, codes, and safety requirements before practical application.

Temperature Profile of Heating Device using Data Logger and Scilab

Topic: Heating Device Temperature Profiling and Analysis;

Subject: Heat Transfer;

Tool: Temperature Data Logger and Scilab;

By: Gani Comia, Dec. 2024;

Data logger nowadays has a built-in software and computer interface to visualize in real-time as in this case the temperature profile of a device being monitored. However, for data analysis and reporting, a software like Scilab comes into picture.

Presented in this article is an example of visual report of temperature data over time. Data points from the logger can be extracted in the form of .CSV file. Transforming this data may sometimes be needed before an analysis be taken. The Scilab scripting facilitates the analysis and presentation of data for reporting.

Fig.1. Temperature Profile using Scilab

Shown below is the script to generate the graph of the data points.

Scilab Script:

// Copyright (C) 2024 - Gani Comia
// Date of creation: 6 Dec 2024
// Heating Oven Temperature Profile using Data Logger & Scilab
clear;clc;
// data extraction
clear importdata;

function [data]=importdata(filename)
    data = csvRead(filename, ",", ".", "double")
endfunction

[A] = importdata("temp_profile_after_repair_xform.csv");

// data transformation
time = A(2:$,13);
ovenWallTemp = A(2:$,2);
partTemp = A(2:$,9);
// average temperature
avgWallTemp = mean(ovenWallTemp)
avgPartTemp = mean(partTemp)

// data visualization
clf;
f=gcf();
f.figure_size=[900,600];
plot(time,ovenWallTemp,"-b","linewidth",1.75)
plot(time,partTemp,"--r","linewidth",1)
legend(["Oven Wall","Workpiece"],"in_upper_right")
title("Heating Oven Temp Profile from Data Logger")
xlabel("Time [ min ]")
ylabel("Temperature [ deg C ]")

a=gca(); 
// x & y-axis range
a.data_bounds = [0 230;780 250]; 

// additional plot properties
t = 0:60:840;
temp = [230 232 250];

// hourly vertical lines
for i = 2:length(t)
    plot([t(i),t(i)],[temp(1),temp(3)],"--y","linewidth",1)
end

// target temperature horizontal lines
plot([t(1),t($)],[temp(2),temp(2)],"--k","linewidth",1)
xstring(5,232,"Std = 232 C")
xstring(625,232,"https://gani-mech-toolbox.blogspot.com")

// hourly label 
minute = 60;
hour = 1;
for i = 2:length(t)
    note = [string(hour),"hour(s)"]
    xstring(minute,235,note,-90)
    minute = minute + 60;
    hour = hour + 1;
end

// average oven temp lines
plot([t(1),t($)],[avgWallTemp,avgWallTemp],"--b","linewidth",1)
noteWall = ["Avg Wall Temp:",string(round(avgWallTemp)),"C"]
xstring(5,avgWallTemp+0.8,noteWall)

// average part temp lines
plot([t(1),t($)],[avgPartTemp,avgPartTemp],"--r","linewidth",1)
notePart = ["Avg Part Temp:",string(round(avgPartTemp)),"C"]
xstring(5,avgPartTemp+0.5,notePart)

This toolbox of Scilab script can be handy reference to analyze data points and present it in the form of graphs for better understanding of numerical data.

Feel free to comment for inquiry, clarification, or suggestion for improvement. Drop your email to request the softcopy of the file.

Disclaimer: The formulas and calculations presented are for technical reference only. Users must verify the accuracy and ensure compliance with applicable engineering standards, codes, and safety requirements before practical application.

CNC Machining Feed Rate for Material Surface Finish in Ra using Scilab

Topic: Turning (Lathe) Feed Rate Calculation for a given Average Roughness;

Subject: Machine Shop Technology;

Tool: Scilab;

By: Gani Comia, Nov. 2024;

Average Roughness, or Ra, is a parameter to measure the surface finish of a workpiece. Ra is typically expressed in micrometers \( \mu m \). Shown here are some of the common Ra values for different surface finishes:

• Polished; Ra 0.05~0.2 \( \mu m \)
• Machined (Finish): Ra 0.8~3.2 \( \mu m \)
• Machined (Rough): Ra 3.2~12.5 \( \mu m \)
• Casted: Ra 12.5~50 \( \mu m \)

In turning metal workpieces, feed rate, f (mm/rev), is one of the considerations to achieve the required surface finish. Let us say we wanted to achieve an Ra 3.2 \( \mu m \) and we have two types of insert, round and rhombic. The machining parameter, which is feed rate, can be calculated and used in the CNC turning program. Below is the basis for calculating feed rate in mm/rev and its Scilab script.

$$f\;=\;{2}\;{\sqrt{\frac{r\;R_a}{125}}}$$

Scilab Script

// Machining Feed Rate Calculation in CNC Turning

// Gani Comia, Jan. 2023

clear;clc;

// given insert type and corner radius
// round corner, r=6.35 mm, rhombic corner, r=0.2 mm
r = [6.35 0.2];                                 // mm
// standard Ra, surface finish
Ra_std = [0.1 0.2 0.4 0.8 1.6 3.2 6.3 12.5];    // um
Ra = linspace(0.1,12.5,200);                    // um

// formula feed/rev for insert radius and Ra requirement
function f=feedPerRev(r, Ra)
    f = 2.*sqrt(r.*Ra./125);  // mm/rev, feed per rev
endfunction

// calculation of feed rate for standard Ra
f_std_round = feedPerRev(r(1), Ra_std);
f_std_rhombic = feedPerRev(r(2), Ra_std);
// table of feed rate for standard Ra
mprintf("\n Ra(um)\t\tf(mm/rev) Round\tf(mm/rev) Rhombic\n")
Table = [Ra_std' f_std_round' f_std_rhombic'];
mprintf("\n %3.2f \t\t %3.3f \t\t %3.3f\n", Table)

// calculation of feed rate from 0.1 to 12.5 um, Ra
f_round = feedPerRev(r(1), Ra);
f_rhombic = feedPerRev(r(2), Ra);

// plotting results
clf;
plot(Ra, f_round, "b-", "linewidth", 1.5)
plot(Ra, f_rhombic, "r-", "linewidth", 1.5)

title("CNC Lathe Feed Rate for Surface Finish")
xlabel(["Average Roughness", "$\Large{R_a\,\;(\mu\,m)}$"])
ylabel(["Feed Rate", "$\Large{f\,\;(mm/rev)}$"])
legend(["Round r = 6.35 mm", "Rhombic r = 0.2 mm"],2)

// case scenario for Ra=3.2 um
Ra_32 = 3.2;
f_round_32 = feedPerRev(r(1), Ra_32);
f_rhombic_32 = feedPerRev(r(2), Ra_32);

// line plot of special concern
xpt_ver = [3.2 3.2]; ypt_ver = [0 f_round_32];
plot(xpt_ver, ypt_ver, "g--")
xpt_hor = [3.2 0]; ypt_hor = [f_round_32 f_round_32];
plot(xpt_hor, ypt_hor, "g--")
xpt_hor_1 = [3.2 0]; ypt_hor_1 = [f_rhombic_32 f_rhombic_32]
plot(xpt_hor_1, ypt_hor_1, "g--")

// plotting points
plot(Ra_32, f_round_32, ".b")
xstring(0.1, f_round_32, ["$\large{f\;=\;0.806\;mm/rev}$"])
plot(Ra_32, f_rhombic_32, ".r")
xstring(0.3, f_rhombic_32, ["$\large{f\;=\;0.143\;mm/rev}$"])
xstring(3.2, 0.30, ["$\Large{R_a\;=\,3.2\;\mu\;m}$"], -90)

Visualization of the relationship of f and Ra for a given nose radius, r, in turning machining is a helpful guide for machinist reference. 

Fig. 1. Turning Feed Rate for Insert Type and Surface Finish.

This toolbox can be your handy reference to calculate the turning feed rate as the initial machining parameter in the CNC program.

Feel free to comment for inquiry, clarification, or suggestion for improvement. Drop your email to request the softcopy of the file.

Disclaimer: The formulas and calculations presented are for technical reference only. Users must verify accuracy and ensure compliance with applicable engineering standards, codes, and safety requirements before practical application.