%%
%程序功能:三次样条插值多项式
%运行方式:运行主程序main3.m
%程序输出:三次样条插值函数图像、原函数图像
%%
clear all; close all; clc;
%% 原函数初始化,以便与三次样条插值多项式进行对比
syms x
f = 1 ./ (1 + x.^2);
f1d = diff(f, x, 1); %diff(f, x, n)表示f对x求n次导数
f2d = diff(f, x, 2);
%% 给出边界条件,选取I型边界条件
x = -5;
S1d_left = subs(f1d);
x = 5;
S1d_right = subs(f1d);
%% 给出插值节点及其相应函数值
xSample = linspace(-5,5,11); %在区间[-5,5]内等距离选取11个节点
x = xSample;
fSample = subs(f);
n = length(x); %读取采样节点的个数
%% 求取系数矩阵A
lamda = 1/2 * ones(1,n-1);
lamda(1) = 1; %令λ0 = 1
miu = 1/2 * ones(1,n-1); %H(j+1)/( H(j-1) + H(j) )
miu(10) = 1; %令μn = 1, n = 10
delta = zeros(1,n);
A = zeros(n,n);
Lam = 1; %令λ0 = 1
Miu = 1; %令μn = 1, n = 10
for i = 1 : n
for j = 1 : n
if j-i == 1
A(i,j) = lamda(Lam);
Lam = Lam + 1;
end
if i-j == 1
A(i,j) = miu(Miu);
Miu = Miu + 1;
end
end
end
%% 根据公式推导求取delta向量
for i = 1 : n
if i == 1
delta(1) = 6 * (( fSample(2) - fSample(1) ) - S1d_left);
elseif i == n
delta(n) = 6 * (S1d_right - ( fSample(n) - fSample(n-1) ));
else
delta(i) = 6 * ( fSample(i+1) + fSample(i-1) - 2*fSample(i)) / 2;
end
A(i,i) = 2;
end
%% 求解基本方程矩阵Am=d,求取m关系向量
m = inv(A) * delta'; %求解Am=d方法:m =inv(A)*b' = A\d' = inv(A'*A)*A'*b'
syms S1 S2 S3 S4 S5 S6 S7 S8 S9 S10
S=[S1 S2 S3 S4 S5 S6 S7 S8 S9 S10];%插值函数分为10段
syms x
for i = 1 : n-1
S(i) = m(i) / 6*(-5+i-x)^3 + m(i+1) / 6*(x+6-i)^3 ...
+ (fSample(i) - m(i)/6) * (-5+i-x) + (fSample(i+1) - m(i+1)/6) * (x+6-i);
end
%% 显示其三次样条插值多项式的分段表达式
disp('Three cubic spline interpolation polynomials : S(x) =')
disp(S)
%% 绘制原函数和三次样条插值多项式的曲线
figure
x = -5 : 0.01 : 5;
plot(x, subs(f), '-k');
hold on;
plot(xSample, fSample,'ob');
hold off;
hold on;
for i = 1 : n-1
x = -6+i : 0.01 : -5+i;
plot(x, subs(S(i)), '-r');
end
hold off;
legend('f = 1/(1+x^2)曲线','采样点','3 Spline插值曲线')
三次样条插值程序,把节点由10个换成20个,完成的私信我我给转10元,急!!thank啊
- 写回答
- 好问题 0 提建议
- 追加酬金
- 关注问题
- 邀请回答
-
2条回答 默认 最新
- 技术专家团-Joel 2021-10-25 19:57关注
你好,如此便可,只要把x尺度缩小即可,因为你的x取的是-5+i,然而x实际是-5+i/2,这里就有个缩放:
%% %程序功能:三次样条插值多项式 %运行方式:运行主程序main3.m %程序输出:三次样条插值函数图像、原函数图像 %% clear all; close all; clc; %% 原函数初始化,以便与三次样条插值多项式进行对比 syms x f = 1 ./ (1 + x.^2); f1d = diff(f, x, 1); %diff(f, x, n)表示f对x求n次导数 f2d = diff(f, x, 2); %% 给出边界条件,选取I型边界条件 x = -5; S1d_left = subs(f1d); x = 5; S1d_right = subs(f1d); %% 给出插值节点及其相应函数值 xSample = linspace(-5,5,21); %在区间[-5,5]内等距离选取21个节点 x = xSample; fSample = subs(f); n = length(x); %读取采样节点的个数 %% 求取系数矩阵A lamda = 1/2 * ones(1,n-1); lamda(1) = 1; %令λ0 = 1 miu = 1/2 * ones(1,n-1); %H(j+1)/( H(j-1) + H(j) ) miu(10) = 1; %令μn = 1, n = 10 delta = zeros(1,n); A = zeros(n,n); Lam = 1; %令λ0 = 1 Miu = 1; %令μn = 1, n = 10 for i = 1 : n for j = 1 : n if j-i == 1 A(i,j) = lamda(Lam); Lam = Lam + 1; end if i-j == 1 A(i,j) = miu(Miu); Miu = Miu + 1; end end end %% 根据公式推导求取delta向量 for i = 1 : n if i == 1 delta(1) = 6 * (( fSample(2) - fSample(1) ) - S1d_left); elseif i == n delta(n) = 6 * (S1d_right - ( fSample(n) - fSample(n-1) )); else delta(i) = 6 * ( fSample(i+1) + fSample(i-1) - 2*fSample(i)) / 2; end A(i,i) = 2; end %% 求解基本方程矩阵Am=d,求取m关系向量 m = inv(A) * delta'; %求解Am=d方法:m =inv(A)*b' = A\d' = inv(A'*A)*A'*b' syms S1 S2 S3 S4 S5 S6 S7 S8 S9 S10 S11 S12 S13 S14 S15 S16 S17 S18 S19 S20 S=[S1 S2 S3 S4 S5 S6 S7 S8 S9 S10 S11 S12 S13 S14 S15 S16 S17 S18 S19 S20];%插值函数分为20段 syms x for i = 1 : n-1 S(i) = m(i) / 6*(-5+i-x)^3 + m(i+1) / 6*(x+6-i)^3 ... + (fSample(i) - m(i)/6) * (-5+i-x) + (fSample(i+1) - m(i+1)/6) * (x+6-i); end %% 显示其三次样条插值多项式的分段表达式 disp('Three cubic spline interpolation polynomials : S(x) =') disp(S) %% 绘制原函数和三次样条插值多项式的曲线 figure x = -5 : 0.01 : 5; plot(x, subs(f), '-k'); hold on; plot(xSample, fSample,'ob'); hold off; hold on; for i = 1 : n-1 q = -6+i/2 : 0.01 : -5+i/2; x = -6+i : 0.01 : -5+i; plot((x-5)/2, subs(S(i)), '-r'); % 增加了点,这里就按照正确的增加数除以2即可 end hold off; legend('f = 1/(1+x^2)曲线','采样点','3 Spline插值曲线')
解决 无用评论 打赏 举报
悬赏问题
- ¥15 一道python难题
- ¥15 用matlab 设计一个不动点迭代法求解非线性方程组的代码
- ¥15 牛顿斯科特系数表表示
- ¥15 arduino 步进电机
- ¥20 程序进入HardFault_Handler
- ¥15 oracle集群安装出bug
- ¥15 关于#python#的问题:自动化测试
- ¥20 问题请教!vue项目关于Nginx配置nonce安全策略的问题
- ¥15 教务系统账号被盗号如何追溯设备
- ¥20 delta降尺度方法,未来数据怎么降尺度