How can I create a symbolic function which fits my data so that I can use the symbolic INT function to double integrate?

I have a dataset where the first column specifies the angle and the second one the radius. I need to know the function that calculates the radius as a function of theta. I want to use this function as an input for the symbolic function INT, like in the example 'f' below. The aim is to find the (numerical) area of a shape parameterized by θ evaluating the definite integral over θ between the limits 0 and 2*pi.

data = importdata('hz.data');

data = fliplr(data);

data(:, 1) = data(:, 1).*(pi/180);

syms r theta f

% example function

f =  2*theta^2+3*theta^3+6*theta -12;

Area = eval(int(int(r,r,0, f), theta, 0, 2*pi));

Answer 1

[MathWorks Support Team]

There are many possibilities to for computing A and the area of the hz shape:

1. You can use POLYFIT to calculate the coefficients of the desired polynomial function:

data = importdata('hz.data');

data = fliplr(data);

data(:, 1) = data(:, 1).*(pi/180);

p3 = polyfit(data(:,1),data(:,2),3);

% p3 = 

%      0.0053   -0.0740    0.2619    2.0549

Then you can use the coefficients to define your symbolic function:

syms r theta f

f =  0.0053*theta^3 - 0.0740*theta^2 + 0.2619*theta + 2.0549; % use values of p

A = eval(int(int(r,r,0, f), theta, 0, 2*pi))

2. You can use numerical approximation theory with DBLQUAD (as you requested a double integral) and INTERP1 to solve this issue:

data = importdata('hz.data');

data = fliplr(data);

data(:, 1) = data(:, 1).*(pi/180);

Ai = dblquad(@(r,theta) r.* (r <= interp1(data(:,1),data(:,2),theta,'spline')), 0, max(data(:,2)), 0, 2*pi)

3. You use POLYFIT with different orders (here as an example with order 3 and 7) instead of INTERP1:

data = importdata('hz.data');

data = fliplr(data);

data(:, 1) = data(:, 1).*(pi/180);

p3 = polyfit(data(:,1),data(:,2),3);

theta = 0:0.01:2*pi;

y3 = p3(1)*theta.^3 + p3(2)*theta.^2 + p3(3)*theta + p3(4);

A3 = dblquad(@(r,theta) r.* (r <= p3(1)*theta.^3 + p3(2)*theta.^2 + p3(3)*theta + p3(4)), 0, max(y3), 0, 2*pi)

p7 = polyfit(data(:,1),data(:,2),7);

y7 = p7(1)*theta.^7 + p7(2)*theta.^6 + p7(3)*theta.^5 + p7(4)*theta.^4 + p7(5)*theta.^3 + p7(6)*theta.^2 + ...

      p7(7)*theta + p7(8);

A7 = dblquad(@(r,theta) r.* (r <= p7(1)*theta.^7 + p7(2)*theta.^6 + p7(3)*theta.^5 + p7(4)*theta.^4 + ...

      p7(5)*theta.^3 + p7(6)*theta.^2 + p7(7)*theta + p7(8)), 0, max(y7), 0, 2*pi)