三次样条差值(MATLAB)

September 5, 2009

function sanciytchzh
% made by louis 2009/06/12
promps={'请输入插值点:'};
result=inputdlg(promps,'插值点x=',1,{'[ 0.2000 0.4000 0.6000 0.8000 1.0000]'});
x =str2num(char(result));
promps={'请输入插值点对应的函数值f(x)='};
result=inputdlg(promps,'插值点处函数值',1,{'[0.9798652 0.9177710 0.8080348 0.6386093 0.3843735]'});
y =str2num(char(result));
promps={'请输入边界条件:'};
result=inputdlg(promps,'第一类边界条件',1,{'[ 0.20271,1.5574]'});
s2=str2num(char(result));
promps={'请输入边界条件:'};
result=inputdlg(promps,'第二类边界条件',1,{'[ 0,0]'});
s1=str2num(char(result));
hold on
grid on
plot(x,y,'o');
a1=diyi(x,y,s2);
a2=ziran(x,y,s1);
for i=1:length(a1)
disp(['第一类边界条件的弯矩M',num2str(i)','为:',num2str(a1(i))]);
end
for i=1:length(a1)
disp(['
自然边界条件的弯矩M',num2str(i)','为:',num2str(a2(i))]);
end
legend('插值点','第一类边界条件','自然边界条件');
function M=diyi(x,y,s)
%made by louis 2009/06/12
%第一类边界条件的三次样条插值
g1=chashang(x,y);
h=x(2:length(x))-x(1:length(x)-1);
for i=1:length(h)-1
u(i)=h(i)/(h(i)+h(i+1));
r(i)=1-u(i);
end
g(1)=6/h(1)*(g1(1,1)-s(1));
g(2:length(x))=6*g1(2,:);
g(length(x))=6/h(length(h)-1)*(s(2)-(y(length(y))-y(length(y)-1))/h(length(h)));
A=2*eye(length(g));
for i=2:length(g)-1
A(i,i+1)=r(i-1);
A(i,i-1)=u(i-1);
end
A(1,2)=1;
A(length(g),length(g)-1)=1;
M=zhuigan(A,g);
h(2:length(h)+1)=h;
n=50;%选取50个插值点
x1=[min(x):(max(x)-min(x))/(n-1):max(x)];
for i=1:n
m=x(1);
for k=2:length(x)
if x1(i)<=m+h(k);
y1(i)=M(k-1)*(x(k)-x1(i))^3/(6*h(k))+M(k)*(x1(i)-x(k-1))^3/(6*h(k))+(y(k-1)-M(k-1)*h(k)^2/6)*(x(k)-x1(i))/h(k)+(y(k)-M(k)*h(k)^2/6)*(x1(i)-x(k-1))/h(k);
break;
end
m=m+h(k);
end
end
plot(x1,y1)
function M=ziran(x,y,s2_1)
% made by louis 2009/06/11
% 满足自然条件的三次样条插值
g=chashang(x,y);
g=6*g(2,:);
g(length(x)-1)=[];%g的最后一个零去掉,因为这个0并不是所算出的二阶差商
h=x(2:length(x))-x(1:length(x)-1);
for i=1:length(h)-1
u(i)=h(i)/(h(i)+h(i+1));
r(i)=1-u(i);
end
g(1)=g(1)-u(1)*s2_1(1);
g(length(g))=g(length(g))-r(length(r))*s2_1(1);
A=2*eye(length(g));
for i=1:length(g)-1
A(i,i+1)=r(i);
A(i+1,i)=u(i+1);
end
M(2:length(x)-1)=zhuigan(A,g);
M(length(x))=0;
h(2:length(h)+1)=h;
n=50;%选取50个插值点
x1=[min(x):(max(x)-min(x))/(n-1):max(x)];
for i=1:n
m=x(1);
for k=2:length(x)
if x1(i)<=m+h(k);
y1(i)=M(k-1)*(x(k)-x1(i))^3/(6*h(k))+M(k)*(x1(i)-x(k-1))^3/(6*h(k))+(y(k-1)-M(k-1)*h(k)^2/6)*(x(k)-x1(i))/h(k)+(y(k)-M(k)*h(k)^2/6)*(x1(i)-x(k-1))/h(k);
break;
end
m=m+h(k);
end
end
plot(x1,y1,'r');
function fc=chashang(x,y)%求差商
% made by louis 2009/06/08
format long
clc
n=length(x);
l=length(y)-1;
for i=1:n-1
for j=1:l
fc(i,j)=(y(j+1)-y(j))/(x(j+i)-x(j));
end
clear y;
y=fc(i,:);
l=l-1;
end
function X=zhuigan(A,B)%追赶法解方程组
% made by louis 2009/06/08
clc
[L,U]=lu(A);
n=length(B);
y(1)=B(1);
for i=2:n
y(i)=B(i)-L(i,i-1)*y(i-1);
end
x(n)=y(n)/U(n,n);
for i=n-1:-1:1
x(i)=(y(i)-U(i,i+1)*x(i+1))/U(i,i);
end
X=x;

由于是自己写的,肯定比较麻烦,但肯定能用!由于代码较长刚特意运行了下,现在把结果页贴出来吧。
第一类边界条件的弯矩M1为:-7.4996
第一类边界条件的弯矩M2为:-0.39631
第一类边界条件的弯矩M3为:1.9385
第一类边界条件的弯矩M4为:-16.3111
第一类边界条件的弯矩M5为:50.5842
自然边界条件的弯矩M1为:0
自然边界条件的弯矩M2为:-1.5018
自然边界条件的弯矩M3为:-1.139
自然边界条件的弯矩M4为:-2.8956
自然边界条件的弯矩M5为:0
三次样条插值结果

Comments