数值分析编程题c语言.docx
《数值分析编程题c语言.docx》由会员分享，可在线阅读，更多相关《数值分析编程题c语言.docx（18页珍藏版）》请在冰点文库上搜索。
数值分析编程题c语言
上机实习题一
一、题目：
已知A与b
12.38412,2.115237,-1.061074,1.112336,-0.113584,0.718719,1.742382,3.067813,-2.031743
2.115237,19.141823,-3.125432,-1.012345,2.189736,1.563849,-0.784165,1.112348,3.123124
-1.061074,-3.125432,15.567914,3.123848,2.031454,1.836742,-1.056781,0.336993,-1.010103
1.112336,-1.012345,3.123848,27.108437,4.101011,-3.741856,2.101023,-0.71828,-0.037585
A=-0.113584,2.189736,2.031454,4.101011,19.897918,0.431637,-3.111223,2.121314,1.784137
0.718719,1.563849,1.836742,-3.741856,0.431637,9.789365,-0.103458,-1.103456,0.238417
1.742382,-0.784165,-1.056781,2.101023,-3.111223,-0.103458,14.7138465,3.123789,-2.213474
3.067813,1.112348,0.336993,-0.71828,2.121314,-1.103456,3.123789,30.719334,4.446782
-2.031743,3.123124,-1.010103,-0.037585,1.784317,0.238417,-2.213474,4.446782,40.00001
b={2.1874369,33.992318,-25.173417,0.84671695,1.784317,-86.612343,1.1101230,4.719345,-5.6784392}
1．用Household变换，把A化为三对角阵B（并打印B）。
2．用超松弛法求解BX=b（取松弛因子ω=1.4,X（0）=0,迭代9次）。
3．用列主元素消去法求解BX=b。
二、解题方法的理论依据:
1、用Householder变换的理论依据
﹝1﹞令A0=A,a（ij）1=a（ij）,已知Ar_1即Ar_1=a（ij）r
﹝2﹞Sr=sqrt（pow（a,2））
﹝3﹞a（r）=Sr*Sr+abs（a（r+1,r））*Sr
﹝4﹞y（r）=A（r_1）*u®/a®
﹝5﹞Kr=（/2）*Ur的转置*Yr/a®
﹝6﹞Qr=Yr-Kr*Ur
﹝7﹞Ar=A（r-1）-（Qr*Ur的转置+Ur*Qr的转置）r=1,2,,……,n-2
2、用超松弛法求解
其基本思想：
在高斯方法已求出x（m）,x（m-1）的基础上，组合新的序列，从而加快收敛速度。
其算式：
其中ω是超松弛因子，当ω>1时，可以加快收敛速度
3、用消去法求解
用追赶消去法求Bx=b的方法：
q1[0]=0,u1[0]=0,
x[9]=u1[9]
三、1．计算程序：
#include"math.h"
#include"stdio.h"
#definege8
voidmain（）
{
intsign（doublex）;
doublea[][9]=
{
{12.38412,2.115237,-1.061074,1.112336,-0.113584,0.718719,1.742382,3.067813,-2.031743},
{2.115237,19.141823,-3.125432,-1.012345,2.189736,1.563849,-0.784165,1.112348,3.123124},
{-1.061074,-3.125432,15.567914,3.123848,2.031454,1.836742,-1.056781,0.336993,-1.010103},
{1.112336,-1.012345,3.123848,27.108437,4.101011,-3.741856,2.101023,-0.71828,-0.037585},
{-0.113584,2.189736,2.031454,4.101011,19.897918,0.431637,-3.111223,2.121314,1.784317},
{0.718719,1.563849,1.836742,-3.741856,0.431637,9.789365,-0.103458,-1.103456,0.238417},
{1.742382,-0.784165,-1.056781,2.101023,-3.111223,-0.103458,14.713846,3.123789,-2.213474},
{3.067813,1.112348,0.336993,-0.71828,2.121314,-1.103456,3.123789,30.719334,4.446782},
{-2.031743,3.123124,-1.010103,-0.037585,1.784317,0.238417,-2.213474,4.446782,40.00001}
};
doublek,h,s,w;
inti,j,n,m,g;
doubleu[9],x1[9],y[9],q[9],b1[9][10],x[9];
doubleb[9]=
{2.1874369,33.992318,-25.173417,0.84671695,1.784317,-86.612343,1.1101230,4.719345,-5.6784392};
for（j=0;j<7;++j）/*Household变换*/
{
s=0.0;
for（i=j+1;i<9;++i）
s=s+a[i][j]*a[i][j];
s=sqrt（s）;
h=（a[j+1][j]>0）?
（s*s+s*a[j+1][j]）:
（s*s-s*a[j+1][j]）;
for（g=0;g<9;++g）
{
if（g<=j）
u[g]=0;
elseif（g==j+1）
u[g]=a[j+1][j]+s*sign（a[j+1][j]）;
elseu[g]=a[g][j];
}
for（m=0;m<9;++m）
{
y[m]=0;
for（n=0;n<9;++n）
y[m]=y[m]+a[m][n]*u[n];
y[m]=y[m]/h;
}
k=0;
for（i=0;i<9;++i）
k=k+u[i]*y[i];
k=0.5*k/h;
for（i=0;i<9;++i）
q[i]=y[i]-k*u[i];
for（n=0;n<9;++n）
for（m=0;m<9;++m）
a[m][n]=a[m][n]-（q[m]*u[n]+u[m]*q[n]）;
}
printf（"Household:
\n"）;
for（i=0;i<9;++i）
for（j=0;j<9;++j）
{
if（j%9==0）
printf（"\n"）;
printf（"%-9.5f",a[i][j]）;
}
printf（"\n"）;
w=1.4;/*超松弛法*/
for（i=0;i<9;i++）
x1[i]=0;
for（i=0;i<9;i++）
for（j=0;j<9;j++）
{
if（i==j）
b1[i][j]=0;
elseb1[i][j]=-a[i][j]/a[i][i];
}
for（i=0;i<9;i++）
b1[i][9]=b[i]/a[i][i];
for（n=0;n<9;n++）
for（i=0;i<9;i++）
{
s=0;
for（j=0;j<9;j++）
s=s+b1[i][j]*x1[j];
s=s+b1[i][9];
x1[i]=x1[i]*（1-w）+w*s;
}
for（i=0;i<9;i++）
{
if（i==5）
printf（"\n"）;
printf（"x%d=%-10.6f",i,x1[i]）;
}
printf（"\n"）;
u[0]=a[0][0];/*以下是消去法*/
y[0]=b[0];
for（i=1;i<9;i++）
{
q[i]=a[i][i-1]/u[i-1];
u[i]=a[i][i]-q[i]*a[i-1][i];
y[i]=b[i]-q[i]*y[i-1];
}
x[ge]=y[ge]/u[ge];
for（i=ge-1;i>=0;i--）
x[i]=（y[i]-a[i][i+1]*x[i+1]）/u[i];
for（i=0;i<9;i++）
{
if（i==5）
printf（"\n"）;
printf（"x%d=%-10.6f",i,x[i]）;
}
}
intsign（doublex）
{
intz;
z=（x>=（1e-6）?
1:
-1）;
return（z）;
}
2．打印结果：
Household:
12.38412-4.893080.000000.000000.000000.000000.000000.000000.00000
-4.8930825.398426.494100.000000.000000.000000.000000.000000.00000
0.000006.4941020.611508.243930.000000.000000.000000.000000.00000
0.000000.000008.2439323.42284-13.880070.000000.000000.00000-0.00000
0.000000.000000.00000-13.8800729.698284.534500.000000.000000.00000
0.000000.000000.000000.000004.5345016.006124.881440.000000.00000
0.000000.000000.000000.000000.000004.8814426.01332-4.50363-0.00000
0.000000.000000.000000.000000.000000.00000-4.5036321.254064.50450
0.000000.000000.00000-0.000000.000000.00000-0.000004.5045014.53412
x0=1.073409x1=2.272579x2=-2.856601
x3=2.292514x4=2.112165x5=-6.422586
x6=1.357802x7=0.634259x8=-0.587042
x0=1.075799x1=2.275744x2=-2.855515
x3=2.293099x4=2.112634x5=-6.423838
x6=1.357923x7=0.634244x8=-0.587266
四、问题讨论：
此程序具有很好的通用性。
在GS方法的基础上，已经求出x的第m解，第m-1基础上，经过重新组合得新的序列，而在此新序列收敛速度加快。
上机实习题二
一、题目：
已知函数值如下表：
x
1
2
3
4
5
f（x）
0
0.6931478
1.0986123
1.3862944
1.6094378
x
6
7
8
9
10
f（x）
1.791795
1.9459101
2.079445
2.1972246
2.3025851
f’（x）
f’
（1）=1
f’（0）=0.1
试用三次样条插值求f（4.563）及f′（4.563）的近似值。
二、解题方法的理论依据：
任意划分的三弯矩插值法以及方程组解法中的三对角阵追赶算法。
应用三次样条插值法能够对函数产生很好的逼近效果。
而追赶算法又具有计算量少、方法简单、算法稳定的特点。
方法应用条件：
适用于求复杂函数在给定区间内某一点的函数值，给出函数f（x）在区间[a,b]中的n个插值点，并且给出函数在区间端点处的值。
三、1．计算程序：
#include"stdio.h"
#include"math.h"
#definen11
#definege11
voidmain（）
{
inti,m;
doublee,s,E,q[12],u[12],y[12],c[12],w[12];
doubleb[12]={2,0,4.15888308,6.5916738,8.3177664,9.6566268,10.750557,11.6754606,
12.47667,13.1833476,13.8155106,14.0155106};
doublea[12][12]={{-2,-4},{0},{0},{0},{0},{0},{0},{0},{0},{0},{0},{0}};
a[n][n-1]=4;
a[n][n]=2;
for（i=1;i<11;i++）
{
a[i][i-1]=1;
a[i][i]=4;
a[i][i+1]=1;
}
u[0]=a[0][0];/*消去法求c[i]*/
y[0]=b[0];
for（i=1;i<12;i++）
{
q[i]=a[i][i-1]/u[i-1];
u[i]=a[i][i]-q[i]*a[i-1][i];
y[i]=b[i]-q[i]*y[i-1];
}
c[ge]=y[ge]/u[ge];
for（i=ge-1;i>=0;i--）
c[i]=（y[i]-a[i][i+1]*c[i+1]）/u[i];
printf（"请输入要插的值:
"）;
scanf（"%lf",&E）;
for（i=0;i<12;i++）
{
e=fabs（E-i）;
if（e>=2）
w[i]=0;
elseif（e<=1）
w[i]=0.5*fabs（e*e*e）-e*e+2.0/3.0;
else
w[i]=（-1.0/6.0）*fabs（e*e*e）+e*e-2*fabs（e）+4.0/3.0;
}
s=0.0;
for（i=0;i<12;i++）
s=s+c[i]*w[i];
printf（"f（%lf）=%lf",E,s）;
printf（"\n"）;
printf（"请输入要求的导数的值:
"）;
scanf（"%d",&m）;
printf（"f’（%d）=%lf\n",m,（c[m+1]-c[m-1]）/2.0）;
}
输出结果：
请输入要插的值:
4.563
f（4.563）=1.517932
请输入要求的导数的值:
4.563
f’（4.563）=0.249350
四、问题讨论：
在给均匀分划的插值函数x赋值时，由于使用for循环，误将x[i]=i+1写成x[i]=i，导致运算错误。
此程序具有一定通用性，对于任意划分的三弯矩插值法，只许改动x[i]即可。
求解方程组Mj时，要用到三对角方程组的追赶法（也称Thomas算法）。
变量较多，应注意区分。
求导时注意正负号。
上机实习题三
一、题目：
用Newton法求方程
X7-28x4+14=0
在（0.1,1.9）中的近似根（初始近似值取为区间端点，迭代6次或误差小于0.00001）。
二、解题方法及理论依据：
Newton迭代法是平方收敛于方程f（x）=0在区间[a,b]上的唯一解α，收敛速度较快，循环次数少。
方法应用条件：
ⅰ）f（a）f（b）<0
ⅱ）f”（x）在区间[a,b]上不变号.
ⅲ）f’（x）≠0
ⅳ）|f（c）|/b-a≤|f’（c）|其中c是a,b中使min[|f’（a）,f’（b）]达到的一个,则对任意时近似值x0€[a,b],Newton迭代过程为
xk+1=φ（xk）=x-f（xk）/f’（xk）,k=1,2,3…
算法：
令
故以1.9为起点
三、1．计算程序：
#include"math.h"
main（）
{
floatx,y,f,f1;
scanf（"%f",&x）;
do{
y=x;
f=pow（y,7）-28*pow（y,4）+14;/*定义f（x）的表达式*/
f1=7*pow（y,6）-112*pow（y,3）;/*定义f’（x）的表达式*/
x=y-f/f1;/*Newton迭代法*/
}
while（fabs（x-y）>=1e-5）;/*控制误差小于0.00001*/
printf（"\nTheresultofthequestionis%f\n",x）;
}
2．打印结果：
请输入端点值：
1.90.1
x=0.8454973.030577
四、问题讨论：
程序较为简单。
它的几何意义为xk+1是f（x）在点xk的切线与x轴交点，故也称为切线法，它是平方收敛的，此处取xk=1.9收敛性较好，要注意判断f′（xk）是否为零。
上机实习题四
一、题目：
用Romberg算法求
（允许误差ε=0.00001）。
二、解题方法及理论依据：
龙贝格（Romberg）方法求数值积分
T1（0）=（b-a）/2*[f（a）+f（b）]
T1（l）=（1/2）*{T1（l-1）+（b-a）/2l-1*∑f[a+（2i-1）*（b-a）/2l]}
Tm+1k-1=[4mTm（k）-Tm（k-1）]/（4m-1）
三、1．计算程序：
#include"math.h"
inta=1,b=3;
doublef（doublex）/*求f（x）=3xx1.4（5x+7）sinx2的值*/
{
doublez;
z=pow（3,x）*pow（x,1.4）*（5*x+7）*sin（x*x）;
return（z）;
}
doubles（intl）/*求T1（l）中的（b-a）/2l-1*∑f[a+（2i-1）*（b-a）/2l]*/
{
externa,b;inti,m;doublez=0.0;
m=pow（2,l-1）;
for（i=1;i<=m;i++）z+=f（a+1.0*（2*i-1）/m）;
z*=1.0*（b-a）/m;
return（z）;
}
main（）
{externa,b;
doubleT[20][20];
intm,n,l=0;
T[1][0]=（b-a）/2.0*（f（a）+f（b））;
do/*龙贝格（Romberg）算法*/
{l++;
T[1][l]=0.5*（T[1][l-1]+s（l））;
n=l-1;
for（m=2;n>=0;m++,n--）/*求解Tl（0）*/
T[m][n]=（pow（4,m-1）*T[m-1][n+1]-T[m-1][n]）/（pow（4,m-1）-1.0）;
}
while（fabs（T[l][0]-T[l-1][0]）>=1e-5）;
printf（"\nT[%d][0]=%f",l,T[l][0]）;
}
2．打印结果：
T[8][0]=440.536017
四、问题讨论：
此程序较繁，计算T1（k）需要复化梯形公式，还要用到Richardson外推法，构造新序列，计算新分点的值时，这些数值个数成倍增加。
应用给出所要求的误差ε，当|Tl（0）-Tl+1（0）|<ε时控制循环。
程序具有广泛的通用性。
上机实习题五
一、题目：
用定步长四阶Runge-Kutta法求解
dy1/dt=1
dy2/dt=y3
dy3/dt=1000-1000y2-100y3
y1（0）=0
y2（0）=0
y3（0）=0
h=0.0005,打印yi（0.025）,yi（0.045）,yi（0.085）,yi（0.1）,（i=1,2,3）
二、解题方法及理论依据：
高阶方程组的Runge-Kutta解法
Yn+1=Yn+（1/6）*（K1+2K2+2K3+K4）
K1=h*F（xn,Yn）
K2=h*F（xn+h/2,Yn+K1/2）
K3=h*F（xn+h/2,Yn+K2/2）
K4=h*F（xn+h,Yn+K3）
适用条件：
使用于那些用普通的积分方法解不了的微分方程组.只要知道函数之间的关系和初值就可以不用解出表达式而直接求解函数在要求点的值。
三、1．计算程序：
#include
main（）
{
doublet,h=0.0005,y1=0,y2=0,y3=0,ky1[5],ky2[5],ky3[5];
for（t=h;t<=0.1001;t+=h）/*Runge-Kutta算法具体过程*/
{ky1[1]=h*1;ky2[1]=h*y3;ky3[1]=h*（1000-1000*y2-100*y3）;
ky1[2]=h*1;ky2[2]=h*（y3+ky3[1]/2）;ky3[2]=h*（1000-1000*（y2+ky2[1]/2）-100*（y3+ky3[1]/2））;
ky1[3]=h*1;ky2[3]=h*（y3+ky3[2]/2）;ky3[3]=h*（1000-1000*（y2+ky2[2]/2）-100*（y3+ky3[2]/2））;
ky1[4]=h*1;ky2[4]=h*（y3+ky3[3]）;ky3[4]=h*（1000-1000*（y2+ky2[3]）-100*（y3+ky3[3]））;
y1+=（ky1[1]+2*ky1[2]+2*ky1[3]+ky1[4]）/6;
y2+=（ky2[1]+2*ky2[2]+2*ky2[3]+ky2[4]）/6;
y3+=（ky3[1]+2*ky3[2]+2*ky3[3]+ky3[4]）/6;
升级会员 