ICode9

精准搜索请尝试: 精确搜索
首页 > 其他分享> 文章详细

2021-10-28

2021-10-28 21:33:03  阅读:173  来源: 互联网

标签:yi 10 xi end i1 28 j1 hi 2021


三次样条插值函数matlab代码实现

xi=[27.7 28 29 30]
yi=[4.1 4.3 4.1 3]
si=[3 -4]
diffq(xi,yi)
s=splin3(xi,yi,si)
simplify(s)
%均差表计算
function p=diffq(xi,yi)
n=length(xi);
f=ones(n,n);
%对差商表第一列赋值
for kk=1:n
    f(kk)=yi(kk);
end

for i1=2:n
    for j1=i1:n
        f(j1,i1)=(f(j1,i1-1)-f(j1-1,i1-1))/(xi(j1)-xi(j1-(i1-1)));
    end
end
disp('均差表如下:');
disp(f);
p=f;
end
%三次样条
function s=splin3(xi,yi,si)
p=diffq(xi,yi)
%xi是插值节点横坐标
%yi是插值节点纵坐标
%si是边界条件
n=length(xi);
a=zeros(n,n);
hi=zeros(n-1,1);
%第一种边界条件
u=zeros(n-1,1);
r=zeros(n-1,1);
d=zeros(n,1);
for i=1:n-1
    hi(i)=xi(i+1)-xi(i);
end
u(n-1)=1;
for ii=1:n-2
    u(ii)=(hi(ii)/(hi(ii)+hi(ii+1)));
end
r(1)=1;
for iii=2:n-1
    r(iii)=(hi(iii)/(hi(iii-1)+hi(iii)));
end
d(1)=(6*(p(2,2)-si(1)))/hi(1);
d(n)=(6*(si(2)-(p(n,2))))/hi(n-1);
for i2=2:n-1
    d(i2)=6*p(i2+1,3);
end
a(n,n)=2;
for k=1:n-1
    for j=1:n-1
        if(k==j)
            a(k,j)=2;
            a(k,j+1)=r(k);
            a(k+1,j)=u(j);
        end
    end
end
m=a\d;
syms x
c=zeros(n,1);
b=sym(c);
for i=1:n-1
    b(i,1)=(m(i)*(xi(i+1)-x)^3)/(6*hi(i))+(m(i+1)*(x-xi(i))^3)/(6*hi(i))+((yi(i)-(m(i)*(hi(i)^2)/6))*(xi(i+1)-x)/hi(i))+((yi(i+1)-(m(i+1)*(hi(i)^2)/6))*(x-xi(i))/hi(i));
end 
s=b;  
end

标签:yi,10,xi,end,i1,28,j1,hi,2021
来源: https://blog.csdn.net/qq_45931366/article/details/121023525

本站声明: 1. iCode9 技术分享网(下文简称本站)提供的所有内容,仅供技术学习、探讨和分享;
2. 关于本站的所有留言、评论、转载及引用,纯属内容发起人的个人观点,与本站观点和立场无关;
3. 关于本站的所有言论和文字,纯属内容发起人的个人观点,与本站观点和立场无关;
4. 本站文章均是网友提供,不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属;如您发现该文章侵犯了您的权益,可联系我们第一时间进行删除;
5. 本站为非盈利性的个人网站,所有内容不会用来进行牟利,也不会利用任何形式的广告来间接获益,纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。

专注分享技术,共同学习,共同进步。侵权联系[81616952@qq.com]

Copyright (C)ICode9.com, All Rights Reserved.

ICode9版权所有