计算机上机报告.docx
《计算机上机报告.docx》由会员分享，可在线阅读，更多相关《计算机上机报告.docx（23页珍藏版）》请在冰点文库上搜索。
计算机上机报告
计算方法上机实验报告
上课时间：
2014-2015学年秋学期，6～14周
一．拉格朗日插值------------------------------------------------------1
二．牛顿插值------------------------------------------------------------3
三．改进欧拉法---------------------------------------------------------5
四．四阶龙格-库塔-----------------------------------------------------7
五．牛顿迭代------------------------------------------------------------9
6．复化Simpson公式------------------------------------------------11
七．Romberg算法------------------------------------------------------14
八.Seidel迭代法------------------------------------------------------17
九.Gauss列主元消去法----------------------------------------------20
一.拉格朗日插值
1.程序代码
#include
voidLagrange（）
{
inti=0;
doublea[10],b[10],L,L1,L2,L3,L4,x;
cout<<"x=";
for（i=0;i<4;i++）
{
cin>>a[i];
}
cout<<"y=";
for（i=0;i<4;i++）
{
cin>>b[i];
}
cout<<"x=";
cin>>x;
L1=（x-a[1]）*（x-a[2]）*（x-a[3]）*b[0]/（a[0]-a[1]）/（a[0]-a[2]）/（a[0]-a[3]）;
L2=（x-a[0]）*（x-a[2]）*（x-a[3]）*b[1]/（a[1]-a[0]）/（a[1]-a[2]）/（a[1]-a[3]）;
L3=（x-a[0]）*（x-a[1]）*（x-a[3]）*b[2]/（a[2]-a[0]）/（a[2]-a[1]）/（a[2]-a[3]）;
L4=（x-a[0]）*（x-a[1]）*（x-a[2]）*b[3]/（a[3]-a[0]）/（a[3]-a[1]）/（a[3]-a[2]）;
L=L1+L2+L3+L4;
cout<<"L="<}
voidmain（）
{
Lagrange（）;
cout<}
2.例子
3.运行结果
二．牛顿插值
1.程序代码
#include
#include
voidmain（）
{
intn,i,j;
doubleA[50][50],*x,*y;
cout<<"请输入插值节点数：
";
cin>>n;
x=newdouble[n];
y=newdouble[n];
cout<<"请输入这"<"<for（i=0;i<=n-1;i++）
cin>>x[i]>>y[i];
doubleK=1,xx,N=0,P;
for（i=0;i<=n-1;i++）
{
A[i][0]=x[i];
A[i][1]=y[i];
}
for（j=2;j<=n;j++）
for（i=1;i<=n-1;i++）
{
A[i][j]=（A[i][j-1]-A[i-1][j-1]）/（A[i][0]-A[i-j+1][0]）;
}
for（i=0;i<=n-1;i++）
cout<<"输出第"<"<cout<<"请输入预求值x=";
cin>>xx;
for（i=0;i{
K*=xx-x[i];
N+=A[i+1][i+2]*K;
P=A[0][1]+N;
}
cout<<"Newton插值结果为：
y="<getch（）;
}
2.例子
3.运行结果
三．改进欧拉法
1.程序代码
#include
#include
doublefun（doublex,doubley）
{
return（-0.9*y/（1+2*x））;
}
voidmain（）
{
doublea,b,*y,h,*x,yp,yc;
intn,k;
cout<<"常微分方程为y'=-0.9*y/（1+2*x）"<cout<<"其中0<=x<=1"<cout<<"初值为y（0）=1"<cout<<"请输入计算区间（a,b）：
";
cin>>a>>b;
cout<<"请输入步长h：
";
cin>>h;
cout<<"请输入计算次数：
";
cin>>n;
y=newdouble[n];
x=newdouble[n];
cout<<"请输入初值y（0）=";
cin>>y[0];
x[0]=a;
for（k=0;k<=n;k++）
{
yp=y[k]+h*fun（x[k],y[k]）;
yc=y[k]+h*fun（x[k]+h,yp）;
y[k+1]=0.5*（yp+yc）;
x[k+1]=x[k]+h;
}
cout<<"迭代结果为：
"<for（k=0;k<=n;k++）
cout<<"x"<getch（）;
}
2.例子
3.运行结果
四．四阶龙格-库塔
1.程序代码
#include
#include
doublefun（doublex,doubley）
{
return（x-y）;
}
voidmain（）
{
doublea,b,*y,h,x,k1,k2,k3,k4;
intn,k;
cout<<"常微分方程为y'=x-y"<cout<<"其中0<=x<=1"<cout<<"初值为y（0）=0"<cout<<"请输入计算区间（a,b）：
";
cin>>a>>b;
cout<<"请输入步长h：
";
cin>>h;
cout<<"请输入计算次数：
";
cin>>n;
y=newdouble[n];
cout<<"请输入初值y（0）：
";
cin>>y[0];
x=a;
cout<<"Runge-Kutta迭代法结果为：
"<cout<<"x0="<for（k=0;k<=n-1;k++）
{
k1=fun（x,y[k]）;
k2=fun（x+h/2,y[k]+k1*h/2）;
k3=fun（x+h/2,y[k]+k2*h/2）;
k4=fun（x+h,y[k]+k3*h）;
y[k+1]=y[k]+（h/6）*（k1+2*（k2+k3）+k4）;
cout<<"x"<x=x+h;
}
getch（）;
}
2.例子
3.运行结果
五．牛顿迭代法
1.程序代码（C++代码）
#include
#include
usingnamespacestd;
doublenewtondiedai（doublea,doubleb,doublec,doubled,doublex）;
intmain（）
{
doublea,b,c,d;
doublex=1.5;
cout<<"请依次输入方程四个系数：
";
cin>>a>>b>>c>>d;
x=newtondiedai（a,b,c,d,x）;
cout<return0;
}
doublenewtondiedai（doublea,doubleb,doublec,doubled,doublex）
{
while（abs（a*x*x*x+b*x*x+c*x+d）>0.000001）
{
x=x-（a*x*x*x+b*x*x+c*x+d）/（3*a*x*x+2*b*x+c）;
}
returnx;
}
2.例子
3.运行结果
六．复化Simpson公式
1.程序代码（C++代码）
#include
#include
doublefunction1（doublex）//被积函数
{
doubles;
s=x/（4+x*x）;
returns;
}
doublefunction2（doublex）//被积函数
{
doubles;
s=sqrt（x）;
returns;
}
doubleReiterationOfSimpson（doublea,doubleb,doublen,doublef（doublex））//复化Simpson公式
{
doubleh,fa,fb,xk,xj;
h=（b-a）/n;
fa=f（a）;
fb=f（b）;
doubles1=0.0;
doubles2=0.0;
for（intk=1;k{
xk=a+k*h;
s1=s1+f（xk）;
}
for（intj=0;j{
xj=a+（j+0.5）*h;
s2=s2+f（xj）;
}
doublesn;//和
sn=h/6*（fa+fb+2*s1+4*s2）;//复化Simpson公式
returnsn;
}
main（）
{doublea,b,Result,n;
cout<<"请输入积分下限:
"<cin>>a;
cout<<"请输入积分上限:
"<cin>>b;
cout<<"请输入分割区间数n:
"<cin>>n;
cout<<"复化Simpson公式计算结果:
";
Result=ReiterationOfSimpson（a,b,n,function1）;
cout<}
2.例子（
）
3.运行结果
七．Romberg算法
1.程序代码（C++代码）
#include
#include
usingnamespacestd;
#definef（x）（4/（1+x*x））
#defineepsilon0.0001
#defineMAXREPT10
doubleRomberg（doubleaa,doublebb）
{intm,n;
doubleh,x;
doubles,q;
doubleep;
double*y=newdouble[MAXREPT];
doublep;
h=bb-aa;
y[0]=h*（f（aa）+f（bb））/2.0;
m=1;
n=1;
ep=epsilon+1.0;
while（（ep>=epsilon）&&（m{p=0.0;
for（inti=0;i{x=aa+（i+0.5）*h;
p=p+f（x）;
}
p=（y[0]+h*p）/2.0;
s=1.0;
for（intk=1;k<=m;k++）
{s=4.0*s;
q=（s*p-y[k-1]）/（s-1.0）;
y[k-1]=p;
p=q;
}
p=fabs（q-y[m-1]）;
m=m+1;
y[m-1]=q;
n=n+n;h=h/2.0;
}
return（q）;
}
intmain（）
{doublea,b;
cout<<"Romberg积分,请输入积分范围a,b:
"<cin>>a>>b;
cout<<"积分结果:
"<system（"pause"）;
return0;
}
2.例子（
）
3.运行结果
八．Seidel迭代法
1.程序代码（C++代码）
#include
#include
#definemax100
#defineEPS1e-6
floata[3][3]={{10,-1,-2},{-1,10,-2},{-1,-1,5}};
floatb[3]={7.2,8.3,4.2};
floatx[3]={0,0,0};
floaty[3];
floatS（intm）
{intn;
floatS=0;
floaty;
for（n=0;n<3;n++）
{if（m==n）{}
else
{S+=a[m][n]*x[n];
}
}y=（b[m]-S）/a[m][m];
returny;
}
voidmain（）
{inti;
intF,T=1,k=1;
do
{F=0;
if（T）
{T=0;}
else
{k++;
for（i=0;i<3;i++）
{x[i]=y[i];
}
}
for（i=0;i<3;i++）
{y[i]=S（i）;
}
for（i=0;i<3;i++）
{if（fabs（x[i]-y[i]）>EPS）
{F=1;
}
}printf（"%d\n",k）;
}while（（（F==1）&&k!
=max））;
printf（"迭代次数:
%d\n",k）;
for（i=0;i<3;i++）
{printf（"x[%d]=%1.6f\n",i,y[i]）;
}
}
2.例子（
）
3.运行结果
九．Gauss列主元消去法
1.程序代码（C++代码）
#include
#include
#include
#definemax_dimension20
intn;
staticfloata[max_dimension][max_dimension];
staticfloatb[max_dimension];
staticfloatx[max_dimension];
voidmain（）
{inti;
intj;
intd;
introw;
floattemp;
floatknown_items;
floatl[max_dimension][max_dimension];
printf（"请输入阶数:
"）;
scanf（"%d",&n）;
printf（"\n"）;
printf（"请输入系数矩阵的值:
"）;
printf（"\n"）;
for（i=0;i{printf（"输入第%d行的值:
",i+1）;
for（j=0;j{
scanf（"%f",&a[i][j]）;
}
printf（"\n"）;
}
printf（"请输入常数项的值:
"）;
for（i=0;iscanf（"%f",&b[i]）;
printf（"\n"）;
for（d=0;d{row=d;
for（i=d+1;i{if（fabs（a[i][d]）>fabs（a[row][d]））
row=i;
}
if（row!
=d）
{for（j=d;j{temp=a[row][j];
a[row][j]=a[d][j];
a[d][j]=temp;
}
temp=b[row];
b[row]=b[d];
b[d]=temp;
}
for（i=d+1;i{l[i][d]=-a[i][d]/a[d][d];
for（j=d;j{a[i][j]=a[i][j]+a[d][j]*l[i][d];
}
b[i]=b[i]+b[d]*l[i][d];
}
}
printf（"\n"）;
for（i=n-1;i>-1;i--）
{known_items=0;
for（j=1;j{known_items=known_items+a[i][i+j]*x[i+j];
}
x[i]=（b[i]-known_items）/a[i][i];
}
printf（"方程组的根为:
\n\n"）;
for（i=0;iprintf（"%.5f",x[i]）;
printf（"\n\n"）;
getch（）;
}
2.例子
3.运行结果
升级会员 