三次参数样条插值方法程序设计

1． 次参数样条插值方法程序设计要求：1) 三学习三次参数样条曲线计算原理；2) 对给定的n个已知数据点Pi(xi,yi) i=1,…,n，确定已知端点条件为抛物端。3) 确定计算参数所满足的线性方程组及边界条件4) 设计求解线性组的方法5) 设计绘制生成曲线的方法6) 编写上述过程的程序7) 程序调试、测试，并利用给定坐标绘出图形编写设计报告。

好2处数组下标溢出

1、g数组只有N个元素，下标可访问区间为[0 , N-1]，下面这个代码相当于访问了g[N]了

for(j=N;j>i;j--)

{

g[j]=(g[j]-g[j-1])/(x[j]-x[j-i-1]);

}

2、g数组最大不能到达N，当i=N的时候，tmpg[i]这里会访问溢出

for(i=1;i<=N;i++)

{

tmp=(u-x[i-1]);

newton=newton+tmpg[i];

}

拉格朗日function y=lagrange(x0,y0,x)

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

for i=1:m

z=x(i);

s=00;

for k=1:n

p=10;

for j=1:n

if j~=k

p=p(z-x0(j))/(x0(k)-x0(j));

end

end

s=py0(k)+s;

end

y(i)=s;

end SOR迭代法的Matlab程序

function [x]=SOR_iterative(A,b)

% 用SOR迭代求解线性方程组,矩阵A是方阵

x0=zeros(1,length(b)); % 赋初值

tol=10^(-2); % 给定误差界

N=1000; % 给定最大迭代次数

[n,n]=size(A); % 确定矩阵A的阶

w=1; % 给定松弛因子

k=1;

% 迭代过程

while k<=N

x(1)=(b(1)-A(1,2:n)x0(2:n)')/A(1,1);

for i=2:n

x(i)=(1-w)x0(i)+w(b(i)-A(i,1:i-1)x(1:i-1)'-A(i,i+1:n)x0(i+1:n)')/A(i,i);

end

if max(abs(x-x0))<=tol

fid = fopen('SOR_iter_resulttxt', 'wt');

fprintf(fid,'\n用SOR迭代求解线性方程组的输出结果\n\n');

fprintf(fid,'迭代次数: %d次\n\n',k);

fprintf(fid,'x的值\n\n');

fprintf(fid, '%128f \n', x);

break;

end

k=k+1;

x0=x;

end

if k==N+1

fid = fopen('SOR_iter_resulttxt', 'wt');

fprintf(fid,'\n用SOR迭代求解线性方程组的输出结果\n\n');

fprintf(fid,'迭代次数: %d次\n\n',k);

fprintf(fid,'超过最大迭代次数，求解失败！');

fclose(fid);

end Matlab中龙格-库塔(Runge-Kutta)方法原理及实现龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法。由于此算法精度高，采取措施对误差进行抑制，所以其实现原理也较复杂。该算法是构建在数学支持的基础之上的。龙格库塔方法的理论基础来源于泰勒公式和使用斜率近似表达微分，它在积分区间多预计算出几个点的斜率，然后进行加权平均，用做下一点的依据，从而构造出了精度更高的数值积分计算方法。如果预先求两个点的斜率就是二阶龙格库塔法，如果预先取四个点就是四阶龙格库塔法。一阶常微分方程可以写作：y'=f(x,y),使用差分概念。

(Yn+1-Yn)/h= f(Xn,Yn)推出（近似等于，极限为Yn'）

Yn+1=Yn+hf(Xn,Yn)

另外根据微分中值定理，存在0<t<1,使得

Yn+1=Yn+hf(Xn+th,Y(Xn+th))

这里K＝f(Xn+th,Y(Xn+th))称为平均斜率，龙格库塔方法就是求得K的一种算法。

利用这样的原理，经过复杂的数学推导（过于繁琐省略），可以得出截断误差为O(h^5)的四阶龙格库塔公式：

K1＝f(Xn,Yn);

K2=f(Xn+h/2,Yn+(h/2)K1);

K3=f(Xn+h/2,Yn+(h/2)K2);

K4=f(Xn+h,Yn+hK3);

Yn+1=Yn+h(K1+2K2+2K3+K4)(1/6);

所以，为了更好更准确地把握时间关系，应自己在理解龙格库塔原理的基础上，编写定步长的龙格库塔函数，经过学习其原理，已经完成了一维的龙格库塔函数。

仔细思考之后，发现其实如果是需要解多个微分方程组，可以想象成多个微分方程并行进行求解，时间，步长都是共同的，首先把预定的初始值给每个微分方程的第一步，然后每走一步，对多个微分方程共同求解。想通之后发现，整个过程其实很直观，只是不停的逼近计算罢了。编写的定步长的龙格库塔计算函数：

function [x,y]=runge_kutta1(ufunc,y0,h,a,b)%参数表顺序依次是微分方程组的函数名称，初始值向量，步长，时间起点，时间终点（参数形式参考了ode45函数）

n=floor((b-a)/h);%求步数

x(1)=a;%时间起点

y(:,1)=y0;%赋初值，可以是向量，但是要注意维数

for ii=1:n

x(ii+1)=x(ii)+h;

k1=ufunc(x(ii),y(:,ii));

k2=ufunc(x(ii)+h/2,y(:,ii)+hk1/2);

k3=ufunc(x(ii)+h/2,y(:,ii)+hk2/2);

k4=ufunc(x(ii)+h,y(:,ii)+hk3);

y(:,ii+1)=y(:,ii)+h(k1+2k2+2k3+k4)/6;

%按照龙格库塔方法进行数值求解

end

调用的子函数以及其调用语句：

function dy=test_fun(x,y)

dy = zeros(3,1);%初始化列向量

dy(1) = y(2) y(3);

dy(2) = -y(1) + y(3);

dy(3) = -051 y(1) y(2);

对该微分方程组用ode45和自编的龙格库塔函数进行比较，调用如下：

[T,F] = ode45(@test_fun,[0 15],[1 1 3]);

subplot(121)

plot(T,F)%Matlab自带的ode45函数效果

title('ode45函数效果')

[T1,F1]=runge_kutta1(@test_fun,[1 1 3],025,0,15);%测试时改变test_fun的函数维数，别忘记改变初始值的维数

subplot(122)

plot(T1,F1)%自编的龙格库塔函数效果

title('自编的 龙格库塔函数')

matlab中使用插值函数

插值函数(the function of interpolation )

interp1

调用函数的格式(Syntax)

yi = interp1(x,Y,xi)

yi = interp1(Y,xi)

yi = interp1(x,Y,xi,method)

yi = interp1(x,Y,xi,method,'extrap')

yi = interp1(x,Y,xi,method,extrapval)

pp = interp1(x,Y,method,'pp')

调用格式说明(Description)

yi = interp1(x,Y,xi) 返回矢量X和Y决定的根据输入的节点xi时对应的y的值矢量Y是矢量X的一个函数映射

如果Y是一个矩阵,那么插值结果是一个对应的矩阵

[===================================================

yi = interp1(x,Y,xi) returns vector yi containing elements corresponding to the elements of xi and determined by interpolation within vectors x and Y The vector x specifies the points at which the data Y is given If Y is a matrix, then the interpolation is performed for each column of Y and yi is length(xi)-by-size(Y,2)

===================================================]

yi = interp1(x,Y,xi,method)插值中可以使用的方法: 插值方法 说明 nearest 临近的两点插值 linear 线性插值（默认） spline 三次样条插值 pchip 分段三次Hermite插值多项式插值 cubic （作用于pchip相同） v5cubic 用matlab5版本中断三次样条插值 [====================================================

yi = interp1(x,Y,xi,method) interpolates using alternative methods:

methodDescription

nearestNearest neighbor interpolation

linearLinear interpolation (default)

splinesplineCubic spline interpolation

pchipPiecewise cubic Hermite interpolation

cubic(Same as 'pchip')

v5cubicCubic interpolation used in MATLAB 5

======================================================]

简单程序示例

>>x=[00 01 0195 03 0401 05];

>>y=[039849 039695 039142 038138 036812 035206];

>>plot(x,y);

>>T=interp1(x,y,25,'linear')　%线性插值

(返回结果T=03862)

>> T=interp1(x,y,25,'nearest')　% 两点插值

(返回结果T=03814)

>>T=interp1(x,y,25,'spline')　 % 三次样条插值

(返回结果T =03867)

>>T=interp1(x,y,25,'cubic')　%三次插值

(返回结果T =03867)

问题补充，因字数限制，挪到这

1拉格朗日插值简介：

对给定的n个插值节点x1,x2,…,xn,及其对应的函数值y1=f(x1), y2=f(x2),…, yn=f(xn);使用拉格朗日插值公式，计算在x点处的对应的函数值f(x);

2一维拉格朗日插值c语言程序：

Int lagrange(x0, y0, n, x, y)

Float xo[], yo[], x;

Int n;

Float y

{

Int i, j;

Float p;

y=0;

If (n>1)

{

For(i=0;i<n;i++)

{

P=1;

For(j=1;j<n;j++)

{

If(i!=J)

P=p(x-x0[j]/x0[i]-x0[j]);

}

y=y+py0[i];

Return(0);

}

Else

Return(-1);

}

3例题。已知函数如下表所示，求x=0472处的函数值：

X 046 047 048 049

Y 0484655 04903745 0502750 0511668

计算这个问题的c语言程序如下：

#minclude stdio

#includeM<nathh>

Main()

{

Float x0[4]={ 046, 047,048,049};

Float y0[4]={ 0484655 ,04903745 ,0502750 ,0511668};

Float x, y;

Int n, rtn;

N=4;

X=0472;

Rth=lagrange(x0,y0,n,x,&y);

If(rtn=0)

{

Prinf(“Y(0472)=:%f\n”,y);

}

Else

{

Prinf(“n must be larger than 1\n”);

}

}

计算结果：Y(0472)=:0495553

4问题补充

我的问题与上面的例子类似，计算三维空间一点(x,y,z)对应的函数值（Vx,Vy,Vz）不同的是自变量(point_coordinatetxt)为三维空间散乱点（不是正方体的顶点），因变量(point_datatxt)为矢量（向量 ）。插值算法比较多，常数法，拉格朗日插值，埃特金插值，三阶样条插值等。最简单的就是常数法，查找离目标点（x,y,z）距离最近的已知自变量(Xi,Yi,Zi)，把该点的函数值赋给目标点做函数值，求高手帮忙写写。

以上就是关于三次参数样条插值方法程序设计全部的内容，包括:三次参数样条插值方法程序设计、c++ 牛顿插值法、如何用MATLAB编写的拉格朗日插值算法的程序、二阶龙格-库塔方法的程序和SOR迭代法的程序等相关内容解答，如果想了解更多相关内容，可以关注我们，你们的支持是我们更新的动力！