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

本文格式为Word版，下载可任意编辑三次样条插值课程实践源代码 数值分析编程实践 姓名：邹建雄 学号：2022226055020 主函数 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 数据 — 3 —。