matlab 三弯矩方程_MATLAB三次样条插值之三弯矩法

MATLAB

三次样条插值之三弯矩法

首先说这个程序并不完善，为了实现通用(

1

，

2

，

…

，

n

)格式解题，以及为调用追赶法程序，

没有针对节点数在三个以下的情况进行分类讨论。希望能有朋友给出更好的方法。

首先，通过函数

sanwanj

得到方程的系数矩阵，即追赶法方程的四个向量参数，接下来调用

追赶法(在

intersanwj

函数中)，得到三次样条分段函数系数因子，然后进行多项式合并得

到分段函数的解析式，程序最后部分通过判断输入值的区间自动选择对应的分段函数并计算改

点的值。附：追赶法程序

chase

%%%%%%%%%%%%%%

function [newv,w,newu,newd]=sanwj(x,y,x0,y0,y1a,y1b)

%

三弯矩样条插值

%

将插值点分两次输入，

x0 y0

单独输入

%

边值条件

a

的二阶导数

y1a

和

b

的二阶导数

y1b

，这里建议将

y1a

和

y1b

换成

y2a

和

y2b

，以便于和三转角代码相区别

n=length(x);m=length(y);

if m~=n

error('x or y

输入有误，再来

');

end

v=ones(n-1,1);u=ones(n-1,1);d=zeros(n-1,1);

w=2*ones(n+1);

h0=x(1)-x0;

h=zeros(n-1,1);

for k=1:n-1

h(k)=x(k+1)-x(k);

end

v(1)=h0/(h0+h(1));

u(1)=1-v(1);

d(1)=6*((y(2)-y(1))/h(1)-(y(1)-y0)/h0)/(h0+h(1));

%

for k=2:n-1

v(k)=h(k-1)/(h(k-1)+h(k));

u(k)=1-v(k);

d(k)=6*((y(k+1)-y(k))/h(k)-(y(k)-y(k-1))/h(k-1))/(h(k-1)+h(k));

end

newv=[v;1];

newu=[1;u];

d0=6*((y(1)-y0)/h0-y1a)/h0;