ex3.m

来自「数值分析, 同济大学教材<<现代数值数学和计算>>数值分」· M 代码 · 共 25 行

x = [75 76 77 78 79 80]; 
y = [0 2.768 2.833 2.903 2.976 3.062 3.153 0]; 
y0 = [2.768 2.833 2.903 2.976 3.062 3.153]; 
n = length(x);
pp = csape(x, y, 'second');
[breaks, coefs, npoly, ncoefs, dim] = unmkpp(pp);
m = 50;
steps = (x(n) - x(1)) / m;
t = [x(1) : steps : x(n)];
for j = 1 : n-1
  for i = 1 : m+1
    if (t(i) >= x(j)) & (t(i) <= x(j+1))
      s(i) = polyval(coefs(j, :), t(i) - x(j));
    end
  end
end
plot(x, y0, '*');
hold on;
plot(t, s, 'b');

%calculate polyder
pdiff = polyder(coefs(3, :));
ydiff = polyval(pdiff, 77.5 - 77)