Convert interpolation function to symbolic

Hello,

I have a Nx2 .dat file which I import to Matlab as a Nx2 matrix. This .dat file is roughly a collection of [x,y] coordinates of a curve.

What I am trying to do is to obtain not only the curve that intepolates those point but also its form in symbolic code.

My code is the following

clc
clear all

%Importing .dat file and converting to a nx2 matrix
f = fopen('GITT_bi1g20.dat');
data = textscan(f, '%f %f');
data = cell2mat(data); % convert to matrix from cell array
fclose(f);

%Comparison between different interpolation functions
%Polynomial
x = data(:,1);
y = data(:,2);
p = polyfit(x,y,9);
y1 = polyval(p,x);

%pchip() - Hermitian
xx = x;
p = pchip(x,y,xx);

figure(1)
hold on
plot(data(:,1),data(:,2),'m')
plot(xx,p,'--black')
hold off

My biggest wish is actually write the pchip() that intepolates the data points as a symbolic function. (So I can do some Calculus stuff with it and other symbolic functions I am using)

Thanks in advance.

Answer 1

[Bjorn Gustavsson]

You already get a piecewise polynomial out of pchip. That is a symbolic function - though not on the form expected by the symbolic functions of matlab so not the format you want. Bot look at the form of your final p and you should be able to figure out what to do next.

HTH

[Luiz Ricardo Almeida]

Hi,

Checking the workspace, I see that my p is a vector full of numbers. How should I look (maybe a command?) to it in order to see its piecewise polynomial?

Thank you very much.

[Bjorn Gustavsson]

That's not what I get. What version of matlab are you using?

There is no problem figuring out what's going on for me running this (and checking the help/documentation

for pchip, polyval and polyfit):

x = sort(randn(5,1));
y = sort(randn(5,1))+randn(5,1);
p = pchip(x,y)
% Thisreturns a struct:
% p = 
%
%  struct with fields:
%
%      form: 'pp'
%    breaks: [-2.9443 -1.1471 -1.0689 -0.8095 0.8884]
%     coefs: [4×4 double]
%    pieces: 4
%     order: 4
%       dim: 1
%
clf
plot(x,y,'.-','markersize',18)
hold on
xi = linspace(min(x),max(x),101);
yi = pchip(x,y,xi);
plot(xi,yi,'r')
plot(x,y,'.-','markersize',18)
y1 = polyval(p.coefs(1,:),xi-p.breaks(1));
y2 = polyval(p.coefs(2,:),xi-p.breaks(2));
y3 = polyval(p.coefs(3,:),xi-p.breaks(3));
y4 = polyval(p.coefs(4,:),xi-p.breaks(4));
ax = axis;                    
plot(xi,y1,'g--'),axis(ax)    
plot(xi,y2,'m--'),axis(ax)
plot(xi,y3,'c--'),axis(ax)
plot(xi,y4,'k--'),axis(ax)

HTH

[Luiz Ricardo Almeida]

Hi, Bjorn,

Thank you for your answer. I figured out that calling pchip as pchip(x,y) and pchip(x,y,z) (where x,y and z are the arguments pchip accepts as inputs) is different. I was using the later, and in that format it yields a vector of numbers. The former would give the expected piecewise funtion.

Anyways, the result is as I expected: too cumbersome. It is almost impossible to transcript by hand the piecewise function pchip gives in a symbolic format to, for exaple, diffrentiate or integrate. Especially when it is divided in like 1000 pieces. I simply gave up. I knew from the beginning that trying to extract inporlating funtions in an "analytic" format is close to madness or a too naive hope.

Feeling that my dreams were shattered and not even close to being fulfilled, I will accept your answer as the best one. Thank you very much.

[Bjorn Gustavsson]

Cheers. I thought that since they are polynomials it would be "reasonably" straight-forward to integrate and differentiate - D(a*x^n) = n*a*x^(n-1) etc. I must admit that I spent no time at all thinking about exactly how the polynomials were represented, so it might get tricky in practice. Then on the other hand polynomials are reasonably well suited for numerical integration so you could wihtout much loss turn to that.