[代码资源] 程序有点问题？？？？？？？帮帮忙

[复制链接]

字体大小: 正常放大

潇湘夜雨123456

2 主题	5 听众	1013 积分

升级 1.3%

TA的每日心情

	开心 2016-12-27 11:45

签到天数: 299 天

[LV.8]以坛为家I

群组: Matlab讨论组

群组: 2011年第一期数学建模

群组: 数学建模培训课堂1

群组: 学术交流A

群组: 西安交大数学建模

电梯直达

1^#

发表于 2012-12-15 11:22 |只看该作者 |倒序浏览

|招呼Ta 关注Ta

function [f,f0] = ThrSample1(x,y,y_1,y_n,x0)
%求已知数据点的第一类三次样条差值多项式及其插值点处的值
%已知数据点的x坐标向量：x
%已知数据点的y坐标向量：y
%左端点的一阶导数：y_1
%右端点的一阶导数：y_n
%插值点的x坐标：x0
%求得的三次样条差值多项式：f
%求得的x0处的插值：f0
syms t
f = 0.0000;
f0 = 0.0000;
if(length(x)==length(y))
n = length(x);
else
disp('x和y的维数不相等！');
return;
end             %维数检查
for i=1:n
if(x(i)<=x0)&&(x(i+1)>=x0)
      ihdex = i;
      break;
end
end             %找到x0所在区间
A = diag(2*ones(1,n)); %求解m的系数矩阵
u = zeros(n-2,1);
lamda = zeros(n-1,1);
c = zeros(n,1);
for i = 2:n-1
u(i-1) = (x(i)-x(x-1))/(x(i+1)-x(i-1));
lamda(i) = (x(i+1)-x(i))/(x(i+1)-x(i-1));
c(i) = 3*lamda(i)*(y(i)-y(i-1))/(x(i)-x(i-1))+3*u(i-1)*(y(i+1)-y(i))/(x(i+1)-x(i));
A(i,i+1) = u(i-1);
A(i,i-1) = lamda(i); %构造系数矩阵及向量c
end
c(1) = 2*y_1;
c(n) = 2*y_n;
m = followup(A,c); %用追赶法求解方程组
h = x(index+1) - x(index);  %x0所在区间长度
f = y(index)*(2*(t-x(index))+h)*(t-x(index+1))^2/h/h/h+...
y(index+1)*(2*(x(index+1)-t)+h)*(t-x(index))^2/h/h/h+...
m(index)*(t-x(index))*(x(index+1)-t)^2/h/h-...
m(index+1)*(x(index+1)-t)*(t-x(index))^2/h/h; %x0所在区间的插值函数
f0 = abs(f,'t',x0)  %x0处的插值

zan

程序