急求二维矩阵三次样条插值的程序,要C或C++语言。望有注释!谢谢! 邮箱[email protected]

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

/下面程序是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

/

yi

=

interp1(x,y,xi,method)

已知样本点坐标x,y,求xi处的函数值yi,插值方法是method

method有以下几种:

'nearest'邻近点插值

'linear'线性插值(默认)

'spline'三次样条函数插值

'cubic'三次函数插值

常用的是'spline'和'cubic'

例子:

x

=

0:10;

y

=

sin(x);

xi

=

0:25:10;

yi

=

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

plot(x,y,'o',xi,yi)

使用Lagrange 插值多项式的方法:

首先把下面的代码复制到M文件中,保存成lagran

function [C,L]=lagran(X,Y)

% input - X is a vector that contains a list of abscissas

% - Y is a vector that contains a list of ordinates

% output - C is a matrix that contains the coefficients of the lagrange interpolatory polynomial

%- L is a matrix that contains the lagrange coefficients polynomial

w=length(X);

n=w-1;

L=zeros(w,w);

for k=1:n+1

V=1;

for j=1:n+1

if k~=j

V=conv(V,poly(X(j)))/(X(k)-X(j));

end

end

L(k,:)=V;

end

C=YL;

然后在命令窗口中输入以下内容:

x=[0 025 05 075 1];

y=[0 03104 06177 07886 1];

lagran(x,y)

ans =

33088 -63851 33164 07599 0

得到的数据就是多项式各项的系数,注意最后一个是常数项,即x^0,

所以表达式为:f=33088x^4-63851x^3+33164x^2 +07599x

求面积就是积分求解

>> f=@(x)33088x^4-63851x^3+33164x^2 +07599x;

>> quad(f,0,1)

ans =

05509

这些点肯定是通过这个多项式的!

出错原因:

函数定义只有一个输出参数,而调用时要求返回两个参数,当然就出错了。

修改:

1、把函数最前面的

function f=language(x,y,x0)

改成

function [f,f0]=language (x,y,x0)

2、另外,最后两句也存在问题:

f0=subs(f,'t',x0);

应为

f0=subs(f,t,x0);

而最后一句则不需要,直接删掉即可。

改后调用实例:

>> [f,f0]=language(x,y,16) 

f = 

-799/3125t(t-1)(t-3/2)(t-2)(t-5/2)(t-3)+561/500t(t-1/2)(t-3/2)(t-2)(t-5/2)(t-3)-133/75t(t-1/2)(t-1)(t-2)(t-5/2)(t-3)+3031/2500t(t-1/2)(t-1)(t-3/2)(t-5/2)(t-3)-399/1250t(t-1/2)(t-1)(t-3/2)(t-2)(t-3)+1411/112500t(t-1/2)(t-1)(t-3/2)(t-2)(t-5/2)

 

f0 =

             099957318144

以上就是关于急求二维矩阵三次样条插值的程序,要C或C++语言。望有注释!谢谢! 邮箱[email protected]全部的内容,包括:急求二维矩阵三次样条插值的程序,要C或C++语言。望有注释!谢谢! 邮箱[email protected]、如何用MATLAB 编写interpl插值函数、求一段matlab插值程序等相关内容解答,如果想了解更多相关内容,可以关注我们,你们的支持是我们更新的动力!

欢迎分享,转载请注明来源:内存溢出

原文地址: https://outofmemory.cn/zz/10216636.html

(0)
打赏 微信扫一扫 微信扫一扫 支付宝扫一扫 支付宝扫一扫
上一篇 2023-05-06
下一篇 2023-05-06

发表评论

登录后才能评论

评论列表(0条)

保存