急求二维矩阵三次样条插值的程序,要C或C++语言。望有注释！谢谢！ 邮箱lsmann9@hotmail.com

/下面程序是C语言程序（标准C）/

/ 计算给定M0,Mn值的三次样条插值多项式 /

/给定离散点(11,04),(12,08),(14,165),(15,18),M0=Mn=0,/

/用M关系式构造三次样条插值多项式S(x),计算S(125)。 /

//

#include <stdioh>

#define Max_N 20

main()

{int i,k,n;

double h[Max_N+1],b[Max_N+1],c[Max_N+1],d[Max_N+1],M[Max_N+1];

double u[Max_N+1],v[Max_N+1],yy[Max_N+1],x[Max_N+1],y[Max_N+1];

double xx,p,q,S;

printf("\nPlease input n value:"); /输入插值点数n/

do

{ scanf("%d",&n);

if(n>Max_N)

printf("\nplease re-input n value:");

}

while(n>Max_N||n<=0);

printf("Input x[i],i=0,%d:\n",n-1);

for(i=0;i<n;i++) scanf("%lf",&x[i]);

printf("Input y[i],i=0,%d:\n",n-1);

for(i=0;i<n;i++) scanf("%lf",&y[i]);

printf("\nInput the M0,Mn value:");

scanf("%lf%lf",&M[0],&M[n]);

printf("\nInput the x value:"); /输入计算三次样条插值函数的x值/

scanf("%lf",&xx);

if((xx>x[n-1])||(xx<x[0]))

{printf("Please input a number between %f and %f\n",x[0],x[n-1]);

return;

}

/计算M关系式中各参数的值/

h[0]=x[1]-x[0];

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

{h[i]=x[i+1]-x[i];

b[i]=h[i]/(h[i]+h[i-1]);

c[i]=1-b[i];

d[i]=6((y[i+1]-y[i])/h[i]-(y[i]-y[i-1])/h[i-1])/(h[i]+h[i-1]);

}

/用追赶法计算Mi,i=1,,n-1/

d[1]-=c[1]M[0];

d[n-1]-=b[n-1]M[n];

b[n-1]=0; c[1]=0; v[0]=0;

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

{u[i]=2-c[i]v[i-1];

v[i]=b[i]/u[i];

yy[i]=(d[i]-c[i]y[i-1])/u[i];

}

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

{M[n-i]=yy[n-i]-v[n-i]M[n-i+1];

}

/计算三次样条插值函数在x处的值/

k=0;

while(xx>=x[k]) k++;

k=k-1;

p=x[k+1]-xx;

q=xx-x[k];

S=(pppM[k]+qqqM[k+1])/(6h[k])+(py[k]+qy[k+1])/h[k]-h[k](pM[k]+qM[k+1])/6;

printf("S(%f)=%f\n",xx,S); /输出/

getch();

}

/----------------------------------- End of file ------------------------------------/

/程序输入输出：

Please input n value:4

Input x[i],i=0,3:

11 12 14 15

Input y[i],i=0,3:

04 08 165 18

Input the M0,Mn value: 0 0

Input the x value:125

S(1250000)=1033171

/

不清楚你的具体问题，给你如下的例子。你可以参考下。

§2 插值问题

在应用领域中，由有限个已知数据点，构造一个解析表达式，由此计算数据点之间的函数值，称之为插值。

实例：海底探测问题

某公司用声纳对海底进行测试，在5×5海里的坐标点上测得海底深度的值，希望通过这些有限的数据了解更多处的海底情况。并绘出较细致的海底曲面图。

一、一元插值

一元插值是对一元数据点(xi,yi)进行插值。

1． 线性插值：由已知数据点连成一条折线，认为相临两个数据点之间的函数值就在这两点之间的连线上。一般来说，数据点数越多，线性插值就越精确。

调用格式：yi=interp1(x,y,xi,’linear’) %线性插值

zi=interp1(x,y,xi,’spline’) %三次样条插值

wi=interp1(x,y,xi,’cubic’) %三次多项式插值

说明：yi、zi、wi为对应xi的不同类型的插值。x、y为已知数据点。

例1：已知数据：

x 0 1 2 3 4 5 6 7 8 9 1

y 3 5 1 14 16 19 6 4 8 15 2

求当xi=025时的yi的值。

程序：

x=0:1:1;

y=[3 5 1 14 16 1 6 4 8 15 2];

yi0=interp1(x,y,0025,'linear')

xi=0:02:1;

yi=interp1(x,y,xi,'linear');

zi=interp1(x,y,xi,'spline');

wi=interp1(x,y,xi,'cubic');

plot(x,y,'o',xi,yi,'r+',xi,zi,'g',xi,wi,'k-')

legend('原始点','线性点','三次样条','三次多项式')

结果：yi0 = 03500

要得到给定的几个点的对应函数值，可用：

xi =[ 02500 03500 04500]

yi=interp1(x,y,xi,'spline')

结果：

yi =12088 15802 13454

(二) 二元插值

二元插值与一元插值的基本思想一致，对原始数据点(x,y,z)构造见世面函数求出插值点数据(xi,yi,zi)。

一、单调节点插值函数，即x,y向量是单调的。

调用格式1：zi=interp2(x,y,z,xi,yi,’linear’)

‘liner’ 是双线性插值 (缺省)

调用格式2：zi=interp2(x,y,z,xi,yi,’nearest’)

’nearest’ 是最近邻域插值

调用格式3：zi=interp2(x,y,z,xi,yi,’spline’)

‘spline’是三次样条插值

说明：这里x和y是两个独立的向量，它们必须是单调的。z是矩阵，是由x和y确定的点上的值。z和x,y之间的关系是z(i,:)=f(x,y(i)) z(:,j)=f(x(j),y) 即：当x变化时，z的第i行与y的第i个元素相关，当y变化时z的第j列与x的第j个元素相关。如果没有对x,y赋值，则默认x=1:n, y=1:m。n和m分别是矩阵z的行数和列数。

例2：已知某处山区地形选点测量坐标数据为：

x=0 05 1 15 2 25 3 35 4 45 5

y=0 05 1 15 2 25 3 35 4 45 5 55 6

海拔高度数据为：

z=89 90 87 85 92 91 96 93 90 87 82

92 96 98 99 95 91 89 86 84 82 84

96 98 95 92 90 88 85 84 83 81 85

80 81 82 89 95 96 93 92 89 86 86

82 85 87 98 99 96 97 88 85 82 83

82 85 89 94 95 93 92 91 86 84 88

88 92 93 94 95 89 87 86 83 81 92

92 96 97 98 96 93 95 84 82 81 84

85 85 81 82 80 80 81 85 90 93 95

84 86 81 98 99 98 97 96 95 84 87

80 81 85 82 83 84 87 90 95 86 88

80 82 81 84 85 86 83 82 81 80 82

87 88 89 98 99 97 96 98 94 92 87

其地貌图为：

对数据插值加密形成地貌图。

程序：

x=0:5:5;

y=0:5:6;

z=[89 90 87 85 92 91 96 93 90 87 82

92 96 98 99 95 91 89 86 84 82 84

96 98 95 92 90 88 85 84 83 81 85

80 81 82 89 95 96 93 92 89 86 86

82 85 87 98 99 96 97 88 85 82 83

82 85 89 94 95 93 92 91 86 84 88

88 92 93 94 95 89 87 86 83 81 92

92 96 97 98 96 93 95 84 82 81 84

85 85 81 82 80 80 81 85 90 93 95

84 86 81 98 99 98 97 96 95 84 87

80 81 85 82 83 84 87 90 95 86 88

80 82 81 84 85 86 83 82 81 80 82

87 88 89 98 99 97 96 98 94 92 87];

mesh(x,y,z) %绘原始数据图

xi=linspace(0,5,50); %加密横坐标数据到50个

yi=linspace(0,6,80); %加密纵坐标数据到60个

[xii,yii]=meshgrid(xi,yi); %生成网格数据

zii=interp2(x,y,z,xii,yii,'cubic'); %插值

mesh(xii,yii,zii) %加密后的地貌图

hold on % 保持图形

[xx,yy]=meshgrid(x,y); %生成网格数据

plot3(xx,yy,z+01,'ob') %原始数据用‘O’绘出

2、二元非等距插值

调用格式：zi=griddata(x,y,z,xi,yi,’指定插值方法’)

插值方法有： linear % 线性插值 (默认)

bilinear % 双线性插值

cubic % 三次插值

bicubic % 双三次插值

nearest % 最近邻域插值

例：用随机数据生成地貌图再进行插值

程序：

x=rand(100,1)4-2;

y=rand(100,1)4-2;

z=xexp(-x^2-y^2);

ti=-2:25:2;

[xi,yi]=meshgrid(ti,ti); % 加密数据

zi=griddata(x,y,z,xi,yi);% 线性插值

mesh(xi,yi,zi)

hold on

plot3(x,y,z,'o')

拉格朗日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软件；

然后，打开写好的函数，复制function后的函数名及参数；

第三，在命令窗口黏贴，将参数替换成需要运算的对象，回车；

最后就会呈现出计算结果。

以上就是关于急求二维矩阵三次样条插值的程序,要C或C++语言。望有注释！谢谢！ 邮箱lsmann9@hotmail.com全部的内容，包括:急求二维矩阵三次样条插值的程序,要C或C++语言。望有注释！谢谢！ 邮箱lsmann9@hotmail.com、matlab插值、如何用MATLAB编写的拉格朗日插值算法的程序、二阶龙格-库塔方法的程序和SOR迭代法的程序等相关内容解答，如果想了解更多相关内容，可以关注我们，你们的支持是我们更新的动力！