常用的插值算法程序.doc

1拉格朗日插值的Matlab实践 Matlab中没有现成的拉格朗日插值函数，必须编写一M文件实现拉格朗日插值。 设n个节点数据以数组x0,y0输入（注意Matlab的数组下标从1开始），m个插值点以数 组x输入，输出数组y为m个插值。 编写一个名为lagrange.m的M文件： functiony=lagrange(x0,y0,x); n=length(x0);m=length(x); fori=1:m z=x(i); s=0.0; fork=1:n p=1.0; forj=1:n ifj=k p=p*(z-x0(j)/(x0(k)-x0(j); end end s=p*y0(k)+s; end y(i)=s; end 作者：lq梦里不知是客2008-4-28 08:28 回复此发言 2回复：拉格朗日插值的Matlab实践 #include #include #include typedefstructdata floatx; floaty; Data;/变量x和函数值y的结构 Datad20;/最多二十组数据 floatf(ints,intt)/牛顿插值法，用以返回插商 if(t=s+1) return(dt.y-ds.y)/(dt.x-ds.x); else return(f(s+1,t)-f(s,t-1)/(dt.x-ds.x); floatNewton(floatx,intcount) intn; while(1) coutn; if(n=count-1)/插值次数不得大于count1次 break; else system(cls); /初始化t，y，yt。 floatt=1.0; floaty=d0.y; floatyt=0.0; /计算y值 for(intj=1;j=n;j+) t=(x-dj-1.x)*t; yt=f(0,j)*t; /coutf(0,j)endl; y=y+yt; returny; floatlagrange(floatx,intcount) floaty=0.0; for(intk=0;kcount;k+)/这儿默认为count1次插值 floatp=1.0;/初始化p for(intj=0;jcount;j+) /计算p的值 if(k=j)continue;/判断是否为同一个数 p=p*(x-dj.x)/(dk.x-dj.x); y=y+p*dk.y;/求和 returny;/返回y的值 voidmain() floatx,y; intcount; while(1) coutcount; if(count=20) break;/检查输入的是否合法 system(cls); /获得各组数据 for(inti=0;icount;i+) cout请输入第i+1di.x; cout请输入第i+1di.y; system(cls); coutx; while(1) intchoice=3; cout请您选择使用哪种插值法计算：endl; cout(0):退出endl; cout(1):Lagrangeendl; cout(2):Newtonendl; coutchoice;/取得用户的选择项 if(choice=2) cout你选择了牛顿插值计算方法，其结果为：; y=Newton(x,count);break;/调用相应的处理函数 if(choice=1) cout你选择了拉格朗日插值计算方法，其结果为：; y=lagrange(x,count);break;/调用相应的处理函数 if(choice=0) break; system(cls); cout输入错误!endl; coutx,yendl;/输出最终结果 作者：lq梦里不知是客2008-4-28 08:34 回复此发言 3回复：拉格朗日插值的Matlab实践 （一）拉格朗日插值多项式 #include #include #include floatlagrange(float*x,float*y,floatxx,intn)/*拉格朗日插值算法*/ inti,j; float*a,yy=0.0;/*a作为临时变量，记录拉格朗日插值多项式*/ a=(float*)malloc(n*sizeof(float); for(i=0;i=n-1;i+) ai=yi; for(j=0;j=20)printf(Error!Thevalueofnmustin(0,20).);getch();return1; if(n=0)printf(Error!Thevalueofnmustin(0,20).);getch();return1; for(i=0;i=n-1;i+) printf(x%d:,i); scanf(%f,&xi); printf(n); for(i=0;i=n-1;i+) printf(y%d:,i);scanf(%f,&yi); printf(n); printf(Inputxx:); scanf(%f,&xx); yy=lagrange(x,y,xx,n); printf(x=%f,y=%fn,xx,yy); getch(); （二）牛顿插值多项式 #include #include #include voiddifference(float*x,float*y,intn) float*f; intk,i; f=(float*)malloc(n*sizeof(float); for(k=1;k=n;k+) f0=yk; for(i=0;i=20)printf(Error!Thevalueofnmustin(0,20).);getch();return1; if(n=0)printf(Error!Thevalueofnmustin(0,20).);getch();return1; for(i=0;i=n-1;i+) printf(x%d:,i); scanf(%f,&xi); printf(n); for(i=0;i=0;i-)yy=yy*(xx-xi)+yi; printf(NewtonInter(%f)=%f,xx,yy); getch(); （三）高斯列主元消去法： #include #include #defineN20 intmain() intn,i,j,k; intmi,tmp,mx; floataNN,bN,xN; printf(nInputn:); scanf(%d,&n); if(nN) printf(Theinputnshouldin(0,N)!n); getch(); return1; if(n=0) printf(Theinputnshouldin(0,N)!n); getch(); return1; printf(Nowinputa(i,j),i,j=0.%d:n,n-1); for(i=0;in;i+) for(j=0;jn;j+) scanf(%f,&aij); printf(Nowinputb(i),i,j=0.%d:n,n-1); for(i=0;in;i+) scanf(%f,&bi); for(i=0;in-2;i+) for(j=i+1,mi=i,mx=fabs(aij);jmx) mi=j; mx=fabs(aji); if(imi) tmp=bi;bi=bmi;bmi=tmp; for(j=i;jn;j+) tmp=aij; aij=amij; amij=tmp; for(j=i+1;jn;j+) tmp=-aji/aii; bj+=bi*tmp; for(k=i;k=0;i-) xi=bi; for(j=i+1;jn;j+) xi-=aij*xj; xi/=aii; for(i=0;in;i+) printf(Answer:nx%d=%fn,i,xi); getch(); return0; #include #include #defineNUMBER20 #defineEsc0x1b #defineEnter0x0d 作者：lq梦里不知是客2008-4-28 08:38 回复此发言 4回复：拉格朗日插值的Matlab实践 floatANUMBERNUMBER+1,ark; intflag,n; exchange(intr,intk); floatmax(intk); message(); main() floatxNUMBER;/*/ intr,k,i,j; charcelect; clrscr(); printf(nnUseGauss.); printf(nn1.JiepleasepressEnter.); printf(nn2.ExitpressEsc.); celect=getch(); if(celect=Esc) exit(0); printf(nninputn=); scanf(%d,&n); printf(nnInputmatrixAandB:); for(i=1;i=n;i+) printf(nnInputa%d1-a%d%dandb%d:,i,i,n,i); /*/ for(j=1;j=n+1;j+)/*/ scanf(%f,&Aij); for(k=1;k=n-1;k+) ark=max(k); if(ark=0)/*/ printf(nnItswrong!);message(); elseif(flag!=k) exchange(flag,k); for(i=k+1;i=n;i+) for(j=k+1;j=1;k-) floatme=0; for(j=k+1;j=n;j+) me=me+Akj*xj; xk=(Akn+1-me)/Akk; for(i=1;i=n;i+) printf(nnx%d=%f,i,xi); message(); exchange(intr,intk)/*/ inti; for(i=1;i=n+1;i+) A0i=Ari; for(i=1;i=n+1;i+) Ari=Aki; for(i=1;i=n+1;i+) Aki=A0i; floatmax(intk)/*/ inti; floattemp=0; for(i=k;itemp) temp=fabs(Aik); flag=i; returntemp; message()/*/ printf(nnGoonEnter,ExitpressEsc!); switch(getch() caseEnter:main(); caseEsc:exit(0); default:printf(nnInputerror!);message(); （四）龙贝格求积公式： #include #include #definef(x)(sin(x)/x) #defineN20 #defineMAX20/*/ #definea2 #defineb4 #definee0.00001/*5*/ 作者：lq梦里不知是客2008-4-28 08:38 回复此发言 5回复：拉格朗日插值的Matlab实践 floatLBG(floatp,floatq,intn) inti; floatsum=0,h=(q-p)/n; for(i=1;in;i+) sum+=f(p+i*h); sum+=(f(p)+f(q)/2; return(h*sum); voidmain() inti; intn=N,m=0; floatTMAX+12; T01=LBG(a,b,n); n*=2; for(m=1;mMAX;m+) for(i=0;im;i+) Ti0=Ti1; T01=LBG(a,b,n); n*=2; for(i=1;i=m;i+) Ti1=Ti-11+(Ti-11-Ti-10)/(pow(2,2*m)-1); if(Tm-11Tm1-e) printf(Answer=%fn,Tm1);getch(); return; （五）牛顿迭代公式： #include #include #include #defineN100 #definePS1e-5 #defineTA1e-5 floatNewton(float(*f)(float),float(*f1)(float),floatx0) floatx1,d=0; intk=0; do x1=x0-f(x0)/f1(x0); if(k+N)|(fabs(f1(x1)PS) printf(nFailed!); getch(); exit(); d=(fabs(x1)PS&fabs(f(x1)TA); returnx1; floatf(floatx) returnx*x*x+x*x-3*x-3; floatf1(floatx) return3.0*x*x+2*x-3; voidmain() floatf(float); floatf1(float); floatx0,y0; printf(Inputx0:); scanf(%f,&x0); printf(x(0)=%fn,x0); y0=Newton(f,f1,x0); printf(nTherootisx=%fn,y0); getch(); （六）牛顿-科特斯求积公式： #include #include intNC(a,h,n,r,f) float(*a); floath; intn,f; float*r; intnn,i; floatds; if(n1000|n2) if(f) printf(nFaild!Checkif1n1000!n,n); return(-1); if(n=2) *r=0.5*(*a)0+(*a)1)*(h); return(0); if(n-4=0) *r=0; *r=*r+0.375*(h)*(*a)n-4+3*(*a)n-3+3*(*a)n-2+(*a)n-1); return(0); if(n/2-(n-1)/2=0) nn=n; else nn=n-3; ds=(*a)0-(*a)nn-1; for(i=2;inn) *r=*r+0.375*(h)*(*a)n-4+3*(*a)n-3+3*(*a)n-2+(*a)n-1); return(0); main() floath,r; intn,ntf,f; inti; floata16; printf(Inputthexi(16):n); for(i=0;i=15;i+) scanf(%d,&ai); h=0.2; f=0; ntf=NC(a,h,n,&r,f); if(ntf=0) printf(nR=%fn,r); else printf(nWrong!Returncode=%dn,ntf); getch(); （七）雅克比迭代： #include #include #defineN20 #defineMAX100 #definee0.00001 intmain() intn; inti,j,k; floatt; floataNN,bNN,cN,gN,xN,hN; printf(nInputdimofn:);scanf(%d,&n); if(nN) printf(Faild!Checkif0nN!n);getch();return1; if(n=0) printf(Faild!Checkif0nN!n);getch();return1; printf(Inputai,j,i,j=0%d:n,n-1); for(i=0;in;i+) for(j=0;jn;j+) scanf(%f,&aij); printf(Inputci,i=0%d:n,n-1); for(i=0;in;i+) scanf(%f,&ci); for(i=0;in;i+) for(j=0;jn;j+) bij=-aij/aii;gi=ci/aii; for(i=0;iMAX;i+) for(j=0;jn;j+) hj=gj; for(k=0;kn;k+) if(j=k)continue;hj+=bjk*xk; t=0; for(j=0;jn;j+) if(tfabs(hj-xj)t=fabs(hj-xj); for(j=0;jn;j+) xj=hj; 作者：lq梦里不知是客2008-4-28 08:38 回复此发言 6回复：拉格朗日插值的Matlab实践 if(te) printf(x_i=n); for(i=0;in;i+) printf(x%d=%fn,i,xi); getch(); return0; printf(after%drepeat,returnn,MAX); getch(); return1; getch(); (八)秦九昭算法： #include floatqin(floata,intn,floatx) floatr=0; inti; for(i=n;i=0;i-) r=r*x+ai; returnr; main() floata50,x,r=0; intn,i; do printf(Inputfrequency:); scanf(%d,&n); while(n1); printf(Inputvalue:); for(i=0;i=n;i+) scanf(%f,&ai); printf(Inputfrequency:); scanf(%f,&x); r=qin(a,n,x); printf(Answer:%f,r); getch(); （九）幂法： #include #include #defineN100 #definee0.00001 #definen3 floatxn=0,0,1; floatann=2,3,2,10,3,4,3,6,1; floatyn; main() inti,j,k; floatxm,oxm; oxm=0; for(k=0;kN;k+) for(j=0;jn;j+) yj=0; for(i=0;in;i+) yj+=aji*xi; xm=0; for(j=0;jxm)xm=fabs(yj); for(j=0;jn;j+) yj/=xm; for(j=0;jn;j+) xj=yj; if(fabs(x