Categories
Programming

Lab 8 – Curve Fitting

Problem 1

In this problem, we are tasked to create a function that either adds or subtracts polynomial vectors, or the coefficients of a polynomial. To do this I simply used a switch statement, and had cases for either “add” or “sub”, as well as accounting for having different sized vectors. To battle that, a temporary vector is made with the larger length then all the elements of the vector are added and the rest is zeros. So to better explain let’s look at the polyadd()function:

%% Polyadd function
function [p] = polyadd(p1,p2,operation)
    if length(p1) == length(p2)
        % If equal then check operation 
        % and do operation
        switch operation 
            case "add" 
                p = p1 + p2;
            case "sub"
                p = p2 - p1;
        end
    else
        % Else add zeros to the shorter one
        if length(p1) < length(p2)
            % Add zeros to p1
            p1_tmp = linspace(0,0,length(p2));
            for i = 1:1:length(p1)
                p1_tmp(i) = p1(i);
            end
            % Have p1 now be p1_tmp
            p1 = p1_tmp;
            % Do operation
            switch operation 
                case "add" 
                    p = p1 + p2;
                case "sub"
                    p = p2 - p1;
            end
        else 
            % Add zeros to p2
            p2_tmp = linspace(0,0,length(p1));
            for i = 1:1:length(p2)
                p2_tmp(i) = p2(i);
            end
            % Have p2 now be p1_tmp
            p2 = p2_tmp;
            % Do operation
            switch operation 
                case "add" 
                    p = p1 + p2;
                case "sub"
                    p = p2 - p1;
            end
        end
    end
end

So let’s test this with different vectors with input:

%% Problem 1
% Get inputs for p1, p2 & operation
p1 = input("Enter p1: \n");
p2 = input("Enter p2: \n");
operation = input("Enter operation (add/sub): \n");
% Get result
p = polyadd(p1,p2,operation);
% Display result
fprintf("Your result is p = ");
disp(p);

%%%%%%%%%%%%%%%%%%%%%%%%%%
Enter p1: 
[9 8 7 6 4]
Enter p2: 
[8 7 6]
Enter operation (add/sub): 
"add"
Your result is p =     17    15    13     6     4
%%%%%%%%%%%%%%%%%%%%%%%%%%%
Enter p1: 
[9 8 7 4 3 2 5 6 8]
Enter p2: 
[0 6 5 4 1 2 33]
Enter operation (add/sub): 
"sub"
Your result is p =     -9    -2    -2     0    -2     0    28    -6    -8

As you can see it works flawlessly, and returns the correct results.

Problem 2

In this problem, we are given the equation \[ x = 2(y-2)^2 + 3 \] & the point \[ P(3,4) \]. With Q being an arbitrary point, and d being the distance between Q & P, we are tasked to do various tasks related to these two variables.

A)

In part A) we are asked to create a formula for the distance, d between P & Q. So on the graph given, we can see Q resides on \[ x= 6 \], so using the formula: \[ dist = \sqrt{(x_2 – x_1)^2 + (y_2 – y_1)^2} \] we can use that to create d:

%% Problem 2
% d(P, Q) = sqrt((x2-x1)^2 + (y2-y1)^2)
% Q = [6, y]
%% A)
P = [3 4];
d = @(y) sqrt((6-P(1)).^2 + (y-P(2)).^2);

B)

In part B) we can use the formula of d to create a d vs y graph in the range of \[ [0,4] \]. So creating a variable y to be the range, all I really had to do was use fplot() with d & y.

The code can be seen below:

%% B) 
y = [0,4];

figure;
fplot(d, y);
title("d versus y graph");
hold off;

We then get the plot below:

Problem 2 Part B) Plot

C)

In part C) we are tasked to find the coordinates of Q when \[ d= 3\] so if the formula for d is: \[ d = \sqrt{(6 – P_x)^2 + (y – P_y) ^2} \] then if we isolate for y, then: \[ y=\pm \sqrt{d^2 – (6-P_x)^2} + P_y \] So with the pretty version seen, let’s see the Matlab version:

%% C)
% Isolate y to calculate values
% Since roots have +- values, 
% y1 = -...
% y2 = +...
% y = sqrt(d^2 - (x2-x1)^2)+4
% Given d = 3
y1 = -1 * sqrt(3^2 - (6-P(1)^2)) + 4;
y2 = sqrt(3^2 - (6-P(1)^2)) + 4;
% Print Coordinates 
fprintf("d is equal to 3 at: \n");
fprintf("(6, %.2f) & (6, %.2f)\n", y1, y2);

Now if we run the code we get the result as:

d is equal to 3 at: 
(6, 0.54) & (6, 7.46)

D)

In part D) we need to find the smallest value of d, so to do so we use a range of 0 to 4 in intervals of 0.25. So if we have a vector, d_vec then we can collect all the values of d and get the smallest using min(). To get the values for d, I created a function dist,so let’s check it out:

%% Make distance formula for D)
function [d] = dist(y)
    d = sqrt((6-3)^2 + (y-4)^2);
end

Now we can see the code for part d):

%% D) 
% To find smallest d, we will have it go through a loop 
% and keep the smallest number
% We will use a range of 0 to 4
d_vec = linspace(0,0, 16);
for i = 1:1:16
     % vector for all d values
    for j = 0:0.25:4
        d_vec(i) = dist(j);
    end
    small_d = min(d_vec); % smallest value
end
% With the smallest value, find y value like c)
d_y1 = -1 * sqrt(small_d^2 - (6-P(1)^2)) + 4;
d_y2 = sqrt(small_d^2 - (6-P(1)^2)) + 4;
% print results
fprintf("The smallest d value is: %d\n", small_d);
fprintf("Q = (6, %.2f) & (6, %.2f)\n", d_y1, d_y2);

Now if we run this code, we get the following result:

The smallest d value is: 3
Q = (6, 0.54) & (6, 7.46)

As we can see, the smallest value of d is 3, which is the same value of d we calculated for Q in Part C).

E)

In part E), we simply need to plot the parabola, and points of Q from part D) and C). First I converted the parabola \[ x = 2(y-2)^2 + 3 \] into \[ y= \pm \sqrt{\cfrac{x-3}{2}} + 2 \] so now we are ready to see the code shown below:

%% E)
% parabola is x = 2(y-2)^2 + 3
% Or: y = +-sqrt((x-3)/2) + 2
parabola1 = @(x) sqrt((x-3)./2) + 2;
parabola2 = @(x) -1 * (sqrt((x-3)./2) + 2);

figure; 
% plot the parabola 1
fplot(parabola1, [0, 12], 'k-');
hold on;
fplot(parabola2, [0, 12], 'k-');
hold on;
% plot points from c)
plot(6,y1,'r*',6,y2,'b*');
hold on;
% plot point from d)
plot(6,d_y1,'m*', 6,d_y2,'g*');
hold off;

If we run the code, we get the plot that shows the parabola and points of Q:

Problem 2 Part E) Plot

Problem 3

In this problem we have a sunflower plant that has data for how long it’s height (in) is for each day, first we are get a cubic curve fit of the data, and estimate the height on day 40, and then use linear & spline interpolation to estimate the height on day 40.

The data we are given are defined in the variables below:

%% Problem 3
Day = [7 21 35 49 63 77 91];
Height =[ 8.5 21 50 77 89 98 99];

A)

To find the cubic curve fit, I used the function polyfit() with the Day and Height variables, and the degree is set as 3 since we are doing a cubic curve fit. With the coefficients, I created an equation and also set the estimated height by “subbing” the values in. With everything set, plotting is simple, and first the graph is plotted, then the curve fit, then the estimated height point. So words can’t say as much as the code does, so let’s see it talk:

%% A) 
% Get curve-fit for 3rd degree 
p = polyfit(Day, Height, 3);
% Make equation with p (x^3, x^2, x, c)
p_eq = @(x) p(1).*x.^3 + p(2).*x.^2 + p(3).*x + p(4);
% To plot for Day 40, we will get the point using p_eq's formula
h_est = p(1) * (40^3) + p(2) * (40^2) + p(3) * 40 + p(4);
figure; 
% Plot graph
plot(Day, Height);
hold on;
% Plot curve-fit
fplot(p_eq, [0, 100]);
hold on;
% Plot Day 40 height estimate 
plot(40, h_est, 'k*');
hold off;

We then get the plot as:

Problem 3 Part A) Plot

We can see that on our graph, we have our curve going close to the graph and we can see our estimated height on Day 40 at around 60, or more accurately 57.56 inches.

B)

In B) we need to find the curve of fit using interpolation methods of spline and linear. With these methods, we will also find the estimated height at Day 40, and to find the interpolation curves I used interp1().

So I first made sample points from x = 7 to 91 in intervals of 0.5, and used that for the spline and linear methods. The estimated height is found having the sample point as x = 40 (Day 40), and with the curve and estimated height, we were ready to graph:

%% B) 
% Find interpolation using spline 
% and linear method
% Let sample point be Day 40 
xi = 7:0.5:91;
% Interpolation curve 
spl = interp1(Day, Height,xi,'spline');
lin = interp1(Day, Height,xi);

% Spline estimate of day 40
spl_est = interp1(Day, Height, 40, 'spline');
lin_est = interp1(Day, Height, 40);

figure;
% Plot graph 
plot(Day, Height, 'k-');
hold on;
% Plot spline iterpolation
plot(xi, spl, 'b--');
hold on;
% Plot linear interpolation
plot(xi, lin, 'm--');
hold on;
% Plot height estimations 
plot(40, spl_est, 'b*');
hold on;
plot(40, lin_est, 'm*');
hold off;

We then get the plot as:

Problem 3 Part B) Plot

As you can see in the plot our points are also close to 60, but there actual values are:

  • Spline: 60.92
  • Linear: 59.64

Which are a lot close or are 60 than in Part A).

Problem 4

In problem 4 we are asked to curve-fit data and use it to estimate the force of the rubber specimen when it’s 11.5 in long. Without furthermore let’s look at part a) and b) for the problem:

A)

In Part A) we are asked to find the curve-fit for the graph using a fourth-order polynomial, and the estimated force when the elongation is 11.5 in. To do this I created a function eq_or_est() to do this for me.

So how this function works is by having three input parameters, the first is the polyfit()vector, its variable name is left as p. Next is the mode parameter, and it determines whether to return the estimate or equation:

  • “est”: returns estimate
  • “eq” : returns equation

The next parameter est is used to calculate what the estimate should be of, in our case we will be doing 11.5. Now we can see how the function looks like:

% Have a function to make the equation and est
function [e] = eq_or_est(p, mode, est)
    switch mode
        case "est" % Finds estimate 
            e = p(1)*est^4 + p(2)*est^3 + p(3)*est^2 + p(4)*est + p(5);
        case "eq" % Makes equation 
e = @(x) p(1).*x.^4 + p(2).*x.^3 + p(3).*x.^2 + p(4).*x + p(5);
    end
end

Now with that function, we are able to find the equation, estimate and plot the results:

%% Problem 4
% Force vector 
force = [0 0.6 0.9 1.16 1.18 1.19 1.24 1.48 1.92 3.12 4.14 5.34 6.22 7.12 7.86 8.42];
% Elongation vector
elongation = [0 1.2 2.4 3.6 4.8 6.0 7.2 8.4 9.6 10.8 12 13.2 14.4 15.6 16.8 18];

%% A)
% Get curve fit for 4th polynomial
p = polyfit(force, elongation, 4);

% Make equation with p, p(x) = ax^4 + bx^3 + cx^2 + dx + e
p_eq = eq_or_est(p, "eq", 0);
% Estimate for 11.5 in
p_est = eq_or_est(p, "est", 11.5);

% Plot for points, curve fit and estimate
figure;
plot(elongation,force,  'k-');
hold on;
fplot(p_eq, [0, 20], 'b--');
hold on;
plot(11.5, p_est, 'm*');
hold off;

When we run this code, we get the plot as:

Problem 4 Part A) Plot

We can see that the estimate can be seen to be in the range from 30 to 50, and it’s actual value is 35.89.

B)

In part B instead of finding the fourth-order curve fit, we need to get the spline interpolation curve fit. To do this, we use the interp1() function again, and we set our sample points from 0 to 18 in 0.5 intervals. With that, we find the spline interpolation for the graph using the sample points, and the estimate with sample point of 11.5. So let’s see in code how that looks like:

%% B) 
% Get spline curve fit
xi = 0:0.5:18;
spline = interp1(elongation, force,xi, 'spline');
% Get spline estimate 
spl_est = interp1(elongation, force, 11.5, 'spline');

% Plot the graph, spline interpolation and estimate
figure;
plot(elongation, force, 'k-');
hold on;
plot(xi, spline, 'b--');
hold on;
plot(11.5, spl_est, 'm*');
hold off;

When we run the code, we get the plot as:

Problem 4 Part B) Plot

As you can see, with the spline interpolation, we get an estimate of 3.72 which is significantly less but most likely more accurate than the fourth-order estimate.

Problem 5

Now I may not be a fan of Python (syntax wise), but it was a better experience than using Matlab, so we are tasked with calculating the GPA’s of an individual with the basis that:

  • A: 4.0
  • B: 3.0
  • C = 2.0
  • D = 1.0
  • F = 0.0

A)

In part A) we need to define the GPA function, so with the GPA array and Hours array we can use a for loop to add up (+ grade * hour) a variable G with an if-elif statement. At first I was tempted to use Python 3.10 match statement, but I felt using what I knew would make things go faster, and I would have to read less documentation. Anyways, once we get this sum, we can divide it by the sum of the hours to get the GPA. So let’s check out how GPA() looks like:

# Problem 5
# A) 
# GPA function
def GPA(Grades, Hours): 
    # define grade values
    # have as float values
    G = 0.0 # will be the sum of the grades * hours
    A = 4.0
    B = 3.0
    C = 2.0
    D = 1.0
    F = 0.0
    for i in range(len(Grades)):
        if Grades[i] == 'A':
            G += A * Hours[i]
        elif Grades[i] == 'B':
            G += B * Hours[i]
        elif Grades[i] == 'C':
            G += C * Hours[i]
        elif Grades[i] == 'D':
            G += D * Hours[i]
        elif Grades[i] == 'F':
            G += F * Hours[i]
    return G / sum(Hours)

B)

Now let’s actually calculate GPA’s, and we are given:

  • grades: B B A D
  • hours: 3 4 3 2

So using that, let’s use GPA() to calculate the GPA using the grades, and hours, and print the results. We can see that done below:

# B) 
grades =  ['B', 'B', 'A','D']
hours = [3, 4,3,2]
# calculate gpas using GPA()
gpa = round(GPA(grades, hours), 2)
print(f'The GPA is : {gpa}')

So when we run python3 lab8_p5.py in a terminal we get:

The GPA is : 2.92
Mustafif

By Mustafif

Admin of website

Leave a Reply

Your email address will not be published.