马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
% exp4_4.m --- II型样条
% 这是我们编的II型样条的程序,供同学参考.
% II型样条是对两个端点附加二阶导数值的样条,也称曲率调整样条
function my_splineII
X=[0 1 2 3];
Y=[0 0.5 2 1.5];
S = spline_II(X,Y); % 调用自编程序,这里两个端点二阶导数为零
fnplt(S); % 作图
hold on
plot(X,Y,'or'); % 画上节点
%----------------- II 型样条 -----------------------
function sp = spline_II(X,Y)
% sp = spline_II(X,Y,Dx2_1,Dx2_n) -- II型三次样条插值
% 两端点二阶导数已知,也称端点曲率调整样条
% 输入: X --- 插值点横坐标(行向量)
% Y --- 插值点纵坐标(行向量)
% Dx2_1 --- 第一个端点的二阶导数值
% Dx2_n --- 最后一个端点的二阶导数值
% 常用Dx2_1=Dx2_n=0即自然样条
% 输出: sp --- 样条(即分段多项式,matlab格式)
n = length(X);
h = diff(X);
d = diff(Y)./h;
d1(2:n-1)=6*diff(d)./(h(1:n-2)+h(2:n-1));
mu(2:n-1)=h(1:n-2)./(h(1:n-2)+h(2:n-1));
la(2:n-1)=1-mu(2:n-1);
%--------------------------------
% 计算三对角方程组
a = mu(3:n-1);
b = 2*ones(n-2,1);
c = la(2:n-2);
u(1:n-2) = d1(2:n-1);
a(n-3)=1-h(n-2)/h(n-1);
b(1)=2+h(1)/h(2);
b(n-2)=2+h(n-2)/h(n-1);
c(1)=1-h(1)/h(2);
a(2:n-2)=a(1:n-3);
a(1)=0;
c(n-2)=0;
%--------------------------------
% 调用追赶法解方程组求得节点上的二阶导数
% 注意这里追赶法与你的下标有所不同
M(2:n-1) = tridiag(a,b,c,u);
M(1)=M(2)*(h(1)+h(2))/h(2)-h(1)/h(2)*M(3);
M(n)=(1+h(n-1)/h(n-2))*M(n-1)-h(n-1)/h(n-2)*M(n-2);
%--------------------------------------
% 下面计算分段多项式的四个系数,以下程序通用,只需模仿即可
% Sk(x) = S(k,1)(x-X(k))^3 + S(k,2)(x-X(k))^2 + S(k,3)(x-X(k)) + S(k,4)
S=zeros(n-1,4);
for k=0:n-2
S(k+1,1)=(M(k+2)-M(k+1))/(6*h(k+1));
S(k+1,2)=M(k+1)/2;
S(k+1,3)=d(k+1)-h(k+1)*(2*M(k+1)+M(k+2))/6;
S(k+1,4)=Y(k+1);
end
sp = mkpp(X,S);
sp.coefs;
% 转换成 Matlab 格式
%-----------------------------------------------
% -------- 求解三对角线性方程组的追赶法 --------
function d = tridiag(a,b,c,d)
n=length(d);
%计算p和q
%b(1)=b(1);
c(1)=c(1)/b(1);
for i = 2:1:n-1
b(i)=b(i)-a(i)*c(i-1);
c(i)=c(i)/b(i);
end
b(n)=b(n)-a(n)*c(n-1);
%求下三角
d(1)=d(1)/b(1);
for i =2:1:n
d(i)=(d(i)-a(i)*d(i-1))/b(i);
end
%求上三角
%d(n)=d(n);
for i =n-1:-1:1
d(i)=d(i)-c(i)*d(i+1);
end
老是错误 Error: Missing variable or function.
怎么搞
谢谢
[ 本帖最后由 lxq 于 2006-10-25 19:58 编辑 ] |