三次样条插值课程实践源代码

数值分析编程实践

姓名：邹建雄 学号：2015226055020

主函数

clear; clc;

X = [0,1,2,3,4,5,6,7,8,9,10];

Y = [2.51,3.30,4.04,4.70,5.22,5.54,5.78,5.40,5.57,5.70,5.80]; dy0 = 0.8;dyn = 0.2;

%X=[1 2 3 4];Y=[1 4 9 16]; %dy0 = 2;dyn = 8;

[h,m] = Cal_mk(X,Y,dy0,dyn);

%Yi=Cal_s(X,Y,h,m,1,1.5); %1.8839 --> 1.8625(自带函数) i=1; d=0.01;

for k = 1:length(X)-1 S(1) = X(1);

for Xi = X(k)+d:d:X(k+1) i = i+1;

S(i) = Cal_s(X,Y,h,m,k,Xi); end end

xi=X(1):d:X(length(X)); plot(X,Y,'o',xi,S); data=[xi;S];

xlswrite('data.xlsx',data);

Cal_mk.m函数

function [h,m] = Cal_mk(X,Y,dy0,dyn) %X为已知数据的横坐标 %Y为已知数据的纵坐标

%mk求得的三次样条插值函数的值 %dy0左端点处的一阶导数 %dyn右端点处的一阶导数

n = length(X)-1; r = zeros(1,n-1); v = zeros(1,n-1); e = zeros(1,n-1); h = zeros(1,n-1); h(1) = X(2)-X(1);

%求解参数获得方程组

for k = 2:n %数组下标从1开始 h(k) = X(k+1)-X(k);

r(k-1) = h(k)/(h(k-1)+h(k)); v(k-1) = h(k-1)/(h(k-1)+h(k));

f0 = (Y(k)-Y(k-1))/h(k-1); f1 = (Y(k+1)-Y(k))/h(k); e(k-1) = 3*(r(k-1)*f0+v(k-1)*f1); end

e(1) = e(1)-r(1)*dy0;

e(n-1) = e(n-1)-v(n-1)*dyn;

%用追赶法求解方程组，解除mk u = zeros(1,n-1); u(1) = 2;

l = zeros(1,n-1); y = zeros(1,n-1); y(1) = e(1); fori = 2:n-1

l(i) = r(i)/u(i-1); u(i) = 2-l(i)*v(i-1); y(i) = e(i)-l(i)*y(i-1); end

x(n-1) = y(n-1)/u(n-1); for t = n-2:-1:1

x(t) = (y(t)-v(t)*x(t+1))/u(t); %m(2)对应x(1)=(y(1)-v(1)*x(2))/u(1) end

m = [dy0 x dyn]; return

Cal_s.m函数

function s = Cal_s(X,Y,h,m,k,x)

s = (x-X(k+1))^2*(h(k) + 2*(x-X(k)))*Y(k)/h(k)^3 + (x-X(k))^2*(h(k) + 2*(X(k+1)-x))*Y(k+1)/h(k)^3 + (x-X(k+1))^2*(x-X(k))*m(k)/h(k)^2 + (x-X(k))^2*(x-X(k+1))*m(k+1)/h(k)^2;

return

数据

本文来源：https://www.bwwdw.com/article/vn1h.html