东南大学计算方法与实习上机实验三

  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
东南大学计算方法与实习
实验报告
学院:电子科学与工程学院
学号:06A12528
姓名:陈毓锋
同组人员:罗关生,肖洋
指导老师:李元庆
1、Lagrange插值多项式的表达式为:
=
但上述公式不适合实际计算:
-对x的求值,每次需要的工作量是O(n2)的;
-增加(或删减)和修改节点,都需要重新计算表达式。
但这只是通常的论断,只要稍作修改,它一样可以适合计算。
#include<cstdlib>
using namespace std;
float CorrctAverageCodeBarycentricLagrange(float x0,float y[],float xx,int n,float h){
int i,k=0;
float lx=1,jc=1,xs=1,*x,*w,*L,*XS,W=0;
cout<<"所用时间为:"<<duration2<<endl;
BL2=BarycentricLagrange(x2,y2,xx2,4);
ACBL=CorrctAverageCodeBarycentricLagrange(0.3,y2,xx2,4,0.1);
cout<<"x="<<setprecision(6)<<xx1<<", L[4]="<<setprecision(7)<<L<<endl;
解析:
1)对Lagrange插值多项式稍作修改:
运用多项式
可以将拉格朗日基本多项式重新写为:
令 ,

则它的优点是当插值点的个数增加一个时,每个ωj都除以(xj-xn+1),就可以得到新的ωn+1,则计算的工作量O(n),比重新计算每个多项式所需要的复杂度O(n2)减少了一个量级。
算法:
1输入xi,yi(i=0,1,2,…,n),令Ln(x)=0;
XS[i]=xs;
}
h=pow(h,n-1);
for(i=0;i<n;i++){
w[i]=L[i]*y[i];
if((n-i)%2==0){
if(i==0){
w[i]*=(-1)/(h*jc);
}
else{
w[i]*=(-1)*XS[i]/(h*jc);
}
}
else if((n-i)%2!=0){
if(i==0){
lx*=xx-x[i];
}
for(i=0;i<n-1;i++){//计算(n-1)!
jc*=n-(i+1);
}
for(i=0;i<n;i++){//将1/(x-xj)看做一个整体L()
L[i]=1/(xx-x[i]);
}
for(i=1;i<n;i++){//(n-1)(n-2)…(n-i)/i
xs*=(float)(n-i)/i;
L[i]=1/(xx-x[i]);
}
for(i=0;i<n;i++){
w[i]=L[i]*y[i];
for(j=0;j<n;j++){
if(j!=i){
w[i]*=1/(x[i]-x[j]);
}
}
W+=w[i];
}
P=lx*W;
return P;
}
对比拉格朗日基本表达式算法情况:
#include<iostream>
cout<<"x="<<setprecision(6)<<xx1<<", BL1[4]="<<setprecision(7)<<BL1<<endl;
cout<<"x="<<setprecision(6)<<xx2<<", BL2[4]="<<setprecision(7)<<BL2<<endl;
cout<<"x="<<setprecision(6)<<xx2<<", ACBL[4]="<<setprecision(7)<<ACBL<<endl;
经过验证,即
此时该均匀节点重新定义的ωj与之前定义的ωj值相同,验证为正确定义。
算法:
3输入xi,yi(i=0,1,2,…,n),令Ln(x)=0;
4对i=1,2,…,n计算
算法程序代码:
#include<iostream>
#include<cmath>
#include<ctime>
#include<iomanip>
float y2[5]={1.2222,1.2681,1.3033,1.3293};
float xx2=0.45,BL2,ACBL;
start1=clock();
do{
a++;
L=Lagrange(x1,y1,xx1,4);
}while(a<=100000);
end1=clock();
duration1=(double)(end1-start1)/CLOCKS_PER_SEC;
在关于稳定性的取值时
在初值+0.0001/rand()取微小量,
可以得出:
则得出该方法的稳定性良好。
通过该程序对该算法的验证,可以看出修正后的拉格朗日插值多项式比拉格朗日基本插值多项式计算量较小,稳定性良好,并且均匀节点时算法也具有这样的优点。
x=new float [n];
w=new float [n];
L=new float [n];
XS=new float [n];
float P=0;
for(i=0;i<n;i++){//给x[i]赋值
x[i]=x0+i*h;
}
for(i=0;i<n;i++){//多项式l(x)=(x-x0)(x-x1)…(x-xn)
float x1[5]={0.4,0.55,0.65,0.8,0.9};//非均匀节点测试
float y1[5]={0.41075,0.57815,0.69675,0.88811,1.02652};
float xx1=0.596,L,BL1;
float x2[4]={0.3,0.4,0.5,0.6};//均匀节点测试
l=new float [n];
for(i=0;i<n;i++){
l[i]=y[i];
for(j=0;j<n;j++){
if(j!=i){
l[i]*=(xx-x[j])/(x[i]-x[j]);
}
}
L+=l[i];
}
delete l;
return L;
}
2)当选择均匀节点时,取xi=x0+ih,可以得到
}
运行结果:
3)实验分析:
利用for循环将拉格朗日基本插值多项式和修正后的拉格朗日插值多项式循环了10万次,通过程序运行如图可以看到分别所用时间为0.22ms和0.248ms,则每个程序所用时间为0.22*10-8s和0.248*10-8s。利用多次计算,得到多次所用时间数据,取平均值
则可以发现修正后的插值多项式比基本多项式运行速度更快一点。另外由于代码段事实上调用进行运算时修正所要经过的循环比基本多,仍能跟基本持平,甚至更快,所以该方法所需要的运算量更快。
w[i]*=1/(h*jc);
}
else{
w[i]*=XS[i]/(h*jc);
}
}
W+=w[i];
}
P=lx*W;
return P;
}
运行结果代码:
void main(){
int a=0,b=0;
clock_t start1,start2,end1,end2;
double duration1,duration2;
2对i=1,2,…,n计算
算法程序代码:
#include<iostream>
#include<cmath>
#include<ctime>
#include<iomanip>
#include<cstdlib>
using namespace std;
float BarycentricLagrange(float x[],float y[],float xx,int n){
#include<cmath>
#include<ctime>
#include<iomanip>
#include<cstdlib>
using namespace std;
float Lagrange(float x[],float y[],float xx,int n){
Biblioteka Baiduint i,j;
float *l,L=0;
int i,j;
float lx=1,*w,*L,W=0;
w=new float [n];
L=new float [n];
float P=0;
for(i=0;i<n;i++){//多项式l(x)=(x-x0)(x-x1)…(x-xn)
lx*=xx-x[i];
}
for(i=0;i<n;i++){//将1/(x-xj)看做一个整体L()
cout<<"所用时间为:"<<duration1<<endl;
start2=clock();
do{
b++;
BL1=BarycentricLagrange(x1,y1,xx1,4);
}while(b<=100000);
end2=clock();
duration2=(double)(end2-start2)/CLOCKS_PER_SEC;
相关文档
最新文档