D03编号241438文件名数值方法程序作业.docx
《D03编号241438文件名数值方法程序作业.docx》由会员分享，可在线阅读，更多相关《D03编号241438文件名数值方法程序作业.docx（32页珍藏版）》请在冰豆网上搜索。
D03编号241438文件名数值方法程序作业
第二章上机作业
.n次拉格朗日插值计算公式Ln（x）
拉格朗日子程序：
doublelaglangre（intn,doublex[],doubley[],doublex1）
{
inti,j;
doubleln=1.0;
doubleLn=0.0;
for（i=0;i<=n;i++）
{
ln=1.0;
for（j=0;j<=n;j++）
{
if（i!
=j）
{
ln=（x1-x[j]）/（x[i]-x[j]）*ln;
}
}
Ln=Ln+ln*y[i];
}
return（Ln）;
}
程序说明：
//总节点数为n+1
//x[]为插值节点
//y[]为对应插值节点的函数值
//x1为求值插值点
计算习题2
（1）
主函数为：
#include
#include
doublemain（）
{intn1=10;
doublexx[11]={1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0};
doubleyy[11]={0,0.0953,0.1823,0.2624,0.3365,0.4055,0.4700,0.5306,0.5877,0.6419,0.6934};
doublexxx=1.54;
doubleyyy=0;
doublelaglangre（intn,doublex[],doubley[],doublex1）;
yyy=laglangre（n1,xx,yy,xxx）;
printf（"Theresultis%1.6f\n",yyy）;
getchar（）;
}
计算结果：
当所求x值为1.54时，
Theresultis0.431806
当所求x值为1.98时，
Theresultis0.683588
总结：
当x值为1.54时，精确值为0.431782，误差为0.000024，当x值为1.98时，精确值为0.683097，误差为0.000491。
2.n次牛顿向前插值计算公式
n次牛顿向前插值计算公式的子程序
doubleNewton（intn,doubleh,doublex[N],doubley[N],doublexx）
{
inti,j;
longints1;
doublet,yy[N][N],s;
doubleNewton=0.0;
doublek=y[0];
t=（xx-x[0]）/h;
s=1.0;
s1=1;
for（i=0;ifor（j=0;j{
yy[i][j]=0.0;
}
for（i=0;ifor（j=0;j{
yy[i][j]=y[j+1]-y[j];
y[j]=yy[i][j];
}
for（i=0;i{
y[i+1]=yy[i][0];
}
y[0]=k;
for（i=0;i<11;i++）
{
if（i==0）
{
s=1.0;
s1=1;
}
else
{
s1=s1*i;
s=s*（t-i+1）;
}
Newton=Newton+s*y[i]/s1;
}
return（Newton）;
}
程序说明：
//n表示插值节点数
//h为步长
//x[N]为插值节点
//y[N]为插值节点对应的插值函数值
//xx为所求插值节点
计算习题2
（1）
主函数为：
#include
#defineN11
doublemain（）
{
doubleNewton（intn,doubleh,doublex[N],doubley[N],doublexx）;
intn1=11;
doublex1[N]={1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0};
doubley1[N]={0,0.0953,0.1823,0.2624,0.3365,0.4055,0.4700,0.5306,0.5877,0.6419,0.6934};
doubleh1=0.1;
doublexx1=1.98;
doubleyy1=Newton（n1,h1,x1,y1,xx1）;
printf（"Theresultis:
%f\n",yy1）;
getchar（）;
}
计算结果：
当所求的x值为1.54时，
Theresultis:
0.431806
当所求的x值为1.98时，
Theresultis:
0.683588
总结：
当x值为1.54时，精确值为0.431782，误差为0.000024，当x值为1.98时，精确值为0.683097，误差为0.000491。
第五章上机作业
．杜利特尔三角分解
杜利特尔三角分解子程序
double*Doolittle（doublea[N][N],doubleb[N]）
{
inti,j,k;
double*p,x[N],y[N],sum;
doublel[N][N],u[N][N];
for（i=0;ifor（j=0;j{
l[i][j]=0;
u[i][j]=0;
}
for（i=0;i{
l[i][i]=1;
u[0][i]=a[0][i];
}
for（i=1;i{
l[i][0]=a[i][0]/u[0][0];
}
for（i=1;i{
for（j=i;j{
sum=0;
for（k=0;k{
sum=sum+l[i][k]*u[k][j];
}
u[i][j]=a[i][j]-sum;
}
for（j=i+1;j{
sum=0;
for（k=0;k{
sum=sum+l[j][k]*u[k][i];
}
l[j][i]=（a[j][i]-sum）/u[i][i];
}
}
for（i=0;i{
sum=0;
if（i==0）
{
y[i]=b[i];
}
else
{
for（k=0;k{
sum=sum+l[i][k]*y[k];
}
y[i]=b[i]-sum;
}
}
printf（"TheresultofLis:
\n"）;
for（i=0;ifor（j=0;j{
printf（"%4.2f",l[i][j]）;
}
printf（"\n"）;
printf（"TheresultofUis:
\n"）;
for（i=0;ifor（j=0;j{
printf（"%4.2f",u[i][j]）;
}
printf（"\n"）;
for（i=N-1;i>=0;i--）
{
sum=0;
if（i==N-1）
{
x[i]=y[i]/u[i][i];
}
else
{
for（k=i;k{
sum=sum+u[i][k]*x[k];
}
x[i]=（y[i]-sum）/u[i][i];
}
}
p=x;
return（p）;
}
程序说明
//该函数设计成返回指针值函数。
//a[N][N]为矩阵A
//b[N]为方程组右端b的值组成的n*1矩阵
//l[N][N]为单位下三角矩阵
//u[N][N]为可逆上三角矩阵
//x[N]为方程组的解，返回指针指向该数组
计算习题五（P179）
//b的初始值设为{1,1,1,1,1,1}
主函数为：
voidmain（）
{
inti,j;
intk=0;
double*q,maxx,c[N];
doublea1[N][N]={1,2,4,7,11,16,2,3,5,8,12,17,4,5,6,9,13,18,7,8,9,10,14,19,11,12,13,14,15,20,16,17,18,19,20,21};
doubleb[N]={1,1,1,1,1,1};
double*Doolittle（doublea[N][N],doubleb[N]）;
doublemax（doublexx[N]）;
for（i=0;i<10;i++）
{
k++;
q=Doolittle（a1,b）;
printf（"k=%d\n",k）;
printf（"Theresultofxis（x[0]...x[N-1]）:
\n"）;
for（j=0;j{
c[j]=*（q+j）;
printf（"%f\n",c[j]）;
}
maxx=max（c）;
printf（"Theresultofbis:
\n"）;
for（j=0;j{
b[j]=c[j]/maxx;
printf（"%f\n",b[j]）;
}
getchar（）;
}
}
其中函数max为：
doublemax（doublex[N]）
{
inti;
doublemax;
max=fabs（x[0]）;
for（i=0;i{
if（max<=fabs（x[i+1]））
max=fabs（x[i+1]）;
}
return（max）;
}
计算结果
TheresultofLis:
1.000.000.000.000.000.002.001.000.000.000.000.004.003.001.000.00
0.000.007.006.001.001.000.000.0011.0010.001.001.001.000.0016.00
15.001.001.001.001.00
TheresultofUis:
1.002.004.007.0011.0016.000.00-1.00-3.00-6.00-10.00-15.000.000.00
-1.00-1.00-1.00-1.000.000.000.00-2.00-2.00-2.000.000.000.000.00
-3.00-3.000.000.000.000.000.00-4.00
k=1
Theresultofxis（x[0]...x[N-1]）:
-1.000000
1.000000
-0.000000
-0.000000
-0.000000
-0.000000
Theresultofbis:
-1.000000
1.000000
-0.000000
-0.000000
-0.000000
-0.000000
k=2
Theresultofxis（x[0]...x[N-1]）:
50.166667
-51.166667
2.000000
0.333333
0.166667
2.500000
Theresultofbis:
0.980456
-1.000000
0.039088
0.006515
0.003257
0.048860
k=3
Theresultofxis（x[0]...x[N-1]）:
-204.964169
108.650380
-3.712269
-0.481542
-0.319218
-2.486971
Theresultofbis:
-1.000000
0.530095
-0.018112
-0.002349
-0.001557
-0.012134
k=4
Theresultofxis（x[0]...x[N-1]）:
468.285487
-160.461507
4.551862
0.409722
0.443818
1.915262
Theresultofbis:
1.000000
-0.342657
0.009720
0.000875
0.000948
0.004090
k=5
Theresultofxis（x[0]...x[N-1]）:
-832.662529
210.162432
-5.170269
-0.185233
-0.554945
-1.679107
Theresultofbis:
-1.000000
0.252398
-0.006209
-0.000222
-0.000666
-0.002017
k=6
Theresultofxis（x[0]...x[N-1]）:
1295.271049
-259.127363
5.872836
-0.164120
0.659122
1.565835
Theresultofbis:
1.000000
-0.200056
0.004534
-0.000127
0.000509
0.001209
k=7
Theresultofxis（x[0]...x[N-1]）:
-1855.320735
308.199805
-6.843905
0.625775
-0.759163
-1.500246
Theresultofbis:
-1.000000
0.166117
-0.003689
0.000337
-0.000409
-0.000809
k=8
Theresultofxis（x[0]...x[N-1]）:
2513.654295
-358.313583
8.224557
-1.192847
0.856489
1.457746
Theresultofbis:
1.000000
-0.142547
0.003272
-0.000475
0.000341
0.000580
k=9
Theresultofxis（x[0]...x[N-1]）:
-3273.096634
410.641679
-10.136369
1.861057
-0.951913
-1.428243
Theresultofbis:
-1.000000
0.125460
-0.003097
0.000569
-0.000291
-0.000436
k=10
Theresultofxis（x[0]...x[N-1]）:
4138.984057
-466.636130
12.690545
-2.627512
1.045951
1.406861
Theresultofbis:
1.000000
-0.112742
0.003066
-0.000635
0.000253
0.000340
2.乔利斯基三角分解方法
乔利斯基三角分解子程序
double*Cholesky（doublea[N][N],doubleb[N]）
{
inti,j,k;
double*p;
doublex[N],y[N];
doublesum,sum2;
doublel[N][N];
doubleu[N][N];
if（a[0][0]<=0）
{
printf（"输入矩阵A不是对称正定矩阵"）;
}
for（i=0;ifor（j=0;j{
l[i][j]=0;
u[i][j]=0;
}
for（i=0;i{
if（i==0）
{
l[i][i]=sqrt（a[i][i]）;
}
l[i][0]=a[i][0]/l[0][0];
}
for（i=1;i{
sum=0;
for（k=0;k{
sum=sum+l[i][k]*l[i][k];
}
l[i][i]=sqrt（a[i][i]-sum）;
for（j=i+1;j{
sum2=0;
for（k=0;k{
sum2=sum2+l[j][k]*l[i][k];
}
l[j][i]=（a[j][i]-sum2）/l[i][i];
}
}
for（i=0;ifor（j=0;j{
u[j][i]=l[i][j];
}
for（i=0;i{
sum=0;
if（i==0）
{
y[i]=b[i]/l[i][i];
}
else
{
for（k=0;k{
sum=sum+l[i][k]*y[k];
}
y[i]=（b[i]-sum）/l[i][i];
}
}
for（i=N-1;i>=0;i--）
{
sum=0;
if（i==N-1）
{
x[i]=y[i]/l[i][i];
}
else
{
for（k=i;k{
sum=sum+u[i][k]*x[k];
}
x[i]=（y[i]-sum）/u[i][i];
}
}
p=x;
return（p）;
}
程序说明：
//N为矩阵的维数
//a[N][N]为方程组A矩阵，为对称正定矩阵
//b[N]为方程组的常数向量
//l[N]为可逆的下三角矩阵
//u[N]为单位上三角矩阵
//x[N]为求解的方程组的解
计算习题9的
（1）（P178）
主函数为：
#include
#include
#defineN4
voidmain（）
{
inti;
double*q,c[N];
doublea1[N][N]={4,2.4,2,3,2.4,5.44,4,5.8,2,4,5.21,7.45,3,5.8,7.45,19.66};
doubleb[N]={12.280,16.928,22.957,50.945};
double*Cholesky（doublea[N][N],doubleb[N]）;
q=Cholesky（a1,b）;
for（i=0;i{
c[i]=*（q+i）;
}
printf（"TheresultofXis:
\n"）;
for（i=0;i{
printf（"%f\n",c[i]）;
}
getchar（）;
}
计算结果：
TheresultofXis:
1.200000
-0.800000
1.700000
2.00000
第六章上机作业
1.雅可比迭代公式
雅可比迭代子程序：
voidJacobi（doublex[N],doublea[N][N],doubleb[N],doubleeps,intCNT）
{
inti,j;
doubley[N],c[N],sum,max;
doubleerr=1.0;
intcount=0;
while（err>eps&&count{
sum=0.0;
count++;
err=0.0;
for（i=0;i{
y[i]=x[i];
}
for（i=0;i{
x[i]=b[i];
for（j=0;j{
if（j!
=i）
x[i]-=a[i][j]*y[j];
}
x[i]/=a[i][i];
c[i]=fabs（x[i]-y[i]）*fabs（x[i]-y[i]）;
sum+=c[i];
}
max=sqrt（sum）;
if（err<=max）
err=max;
}
printf（"TheresultofXis:
\n"）;
for（i=0;i{
printf（"%f\n",x[i]）;
}
printf（"Thenumberofiterationsis:
%d",count）;
}
程序说明：
//x[N]为初值
//a[N][N]为矩阵A，即方程系数矩阵
//b[N]为N*1的b矩阵，即方程右边b的值
//eps为允许误差大小
//CNT为设定的迭代最大次数
//N为矩阵的维数
计算习题15的方程组：
主函数为：
#include
#include
#defineN6
voidmain（）
{
doublex1[N]={1,1,1,1,1,1};
doublea1[N][N]={4,-1,0,-1,0,0,-1,4,-1,0,-1,0,0,-1,4,0,0,-1,-1,0,0,4,-1,0,0,-1,0,-1,4,-1,0,0,-1,0,-1,4};
doubleb1[N]={0,5,0,6,-2,6};
doubleeps1=1.0e-5;
intCNT1=100;
voidJacobi（doublex[N],doublea[N][N],doubleb[N],doubleeps,intCNT）;
Jacobi（x1,a1,b1,eps1,CNT1）;
getchar（）;
}
运行结果为：
TheresultofXis:
1.000000
1.999998
1.000000
1.999998
1.000000
1.999998
Thenumberofiterationsis:
26
2.高斯-塞德尔迭代公式
高斯-塞德尔迭代子程序：
VoidGauss_Seidel（doublex[N],doublea[N][N],doubleb[N],doubleeps,intCNT）
{
inti,j;
doubley[N],c[N],sum,max;
doubleerr=1.0;
intcount=0.0;
while（err>eps&&cou
升级会员 