《数值分析》上机实验报告.docx
《《数值分析》上机实验报告.docx》由会员分享，可在线阅读，更多相关《《数值分析》上机实验报告.docx（22页珍藏版）》请在冰豆网上搜索。
《数值分析》上机实验报告
数值分析上机实验报告
《数值分析》上机实验报告
1.用Newton法求方程
X7-X4+14=0
在（0.1，1.9）中的近似根（初始近似值取为区间端点，迭代6次或误差小于0.00001）。
1.1理论依据：
设函数在有限区间[a，b]上二阶导数存在，且满足条件
令
故以1.9为起点
如此一次一次的迭代，逼近x的真实根。
当前后两个的差<=ε时，就认为求出了近似的根。
本程序用Newton法求代数方程（最高次数不大于10）在（a,b）区间的根。
1.2C语言程序原代码：
#include
#include
main（）
{doublex2,f,f1;
doublex1=1.9;//取初值为1.9
do
{x2=x1;
f=pow（x2,7）-28*pow（x2,4）+14;
f1=7*pow（x2,6）-4*28*pow（x2,3）;
x1=x2-f/f1;}
while（fabs（x1-x2）>=0.00001||x1<0.1）;//限制循环次数
printf（"计算结果：
x=%f\n",x1）;}
1.3运行结果：
1.4MATLAB上机程序
functiony=Newton（f,df,x0,eps,M）
d=0;
fork=1:
M
iffeval（df,x0）==0
d=2;break
else
x1=x0-feval（f,x0）/feval（df,x0）;
end
e=abs（x1-x0）;
x0=x1;
ife<=eps&&abs（feval（f,x1））<=eps
d=1;break
end
end
ifd==1
y=x1;
elseifd==0
y='迭代M次失败';
else
y='奇异'
end
functiony=df（x）
y=7*x^6-28*4*x^3;
End
functiony=f（x）
y=x^7-28*x^4+14;
End
>>x0=1.9;
>>eps=0.00001;
>>M=100;
>>x=Newton（'f','df',x0,eps,M）;
>>vpa（x,7）
1.5问题讨论：
1.使用此方法求方解，用误差来控制循环迭代次数，可以在误差允许的范围内得到比较理想的计算结果。
此程序的不足之处是，所要求解的方程必须满足上述定理的四个条件，但是第二和第四个条件在计算机上比较难以实现。
2.Newton迭代法是一个二阶收敛迭代式，他的几何意义Xi+1是Xi的切线与x轴的交点,故也称为切线法。
它是平方收敛的，但它是局部收敛的，即要求初始值与方程的根充分接近，所以在计算过程中需要先确定初始值。
3.本题在理论依据部分，讨论了区间（0.1,1.9）两端点是否能作为Newton迭代的初值，结果发现0.1不满足条件，而1.9满足，能作为初值。
另外，该程序简单，只有一个循环，且为顺序结构，故采用do-while循环。
当然也可以选择for和while循环。
2.已知函数值如下表：
x
1
2
3
4
5
f（x）
0
0.69314718
1.0986123
1.3862944
1.6094378
x
6
7
8
9
10
f（x）
1.7917595
1.9459101
2.079445
2.1972246
2.3025851
f’（x）
f’
（1）=1
f’（10）=0.1
试用三次样条插值求f（4.563）及f’（4.563）的近似值。
2.1理论依据
这里
，所以只要求出
，就能得出插值函数S（x）。
求
的方法为：
这里
最终归结为求解一个三对角阵的解。
用追赶法解三对角阵的方法如下：
，
综上可得求解方程Ax=d的算法：
2.2C语言程序代码：
#include
#include
voidmain（）
{inti,j,m,n,k,p;
doubleq10,p10,s4,g4,x0,x1,g0=1,g9=0.1;;
doubles[10][10];
doublea[10],b[10],c[10],d[10],e[10],x[10],h[9],u[9],r[9];
doublef[10]={0,0.69314718,1.0986123,1.3862944,1.6094378,
1.7917595,1.9459101,2.079445,2.1972246,2.3025851};
printf（"请依次输入xi:
\n"）;
for（i=0;i<=9;i++）
scanf（"%lf",&e[i]）;//求h矩阵
for（n=0;n<=8;n++）
h[n]=e[n+1]-e[n];
d[0]=6*（（f[1]-f[0]）/h[0]-g0）/h[0];
d[9]=6*（g9-（f[9]-f[8]）/h[8]）/h[8];
for（j=0;j<=7;j++）
d[j+1]=6*（（f[j+2]-f[j+1]）/h[j+1]-（f[j+1]-f[j]）/h[j]）/（h[j]+h[j+1]）;
for（m=1;m<=8;m++）
u[m]=h[m-1]/（h[m-1]+h[m]）;
for（k=1;k<=8;k++）
r[k]=h[k]/（h[k-1]+h[k]）;
for（i=0;i<=9;i++）//求u矩阵
for（p=0;p<=9;p++）
{s[i][p]=0;
if（i==p）s[i][p]=2;}
s[0][1]=1;
s[9][8]=1;
for（i=1;i<=8;i++）
{s[i][i-1]=u[i];
s[i][i+1]=r[i];}
printf（"三对角矩阵为:
\n"）;
for（i=0;i<=9;i++）
for（p=0;p<=9;p++）//求r矩阵
{printf（"%5.2lf",s[i][p]）;
if（p==9）
{printf（"\n"）;}
}
printf（"根据追赶法解三对角矩阵得：
\n"）;
a[0]=s[0][0];
b[0]=d[0];
for（i=1;i<9;i++）
{c[i]=s[i][i-1]/a[i-1];//求d矩阵
a[i]=s[i][i]-s[i-1][i]*c[i];
b[i]=d[i]-c[i]*b[i-1];
if（i==8）
{p10=b[i];
q10=a[i];}}
x[9]=p10/q10;
printf（"M[10]=%lf\n",x[9]）;
for（i=9;i>=1;i--）
{x[i-1]=（b[i-1]-s[i-1][i]*x[i]）/a[i-1];
printf（"M[%d]=%lf\n",i,x[i-1]）;}
printf（"可得s（x）在区间[4,5]上的表达式;\n"）;
printf（"将x=4.563代入得：
\n"）;
x0=5-4.563;
x1=4.563-4;
s4=x[3]*pow（x0,3）/6+x[4]*pow（x1,3）/6+（f[3]-x[3]/6）*（5-4.563）+（f[4]-x[4]/6）*（4.563-4）;
g4=-x[3]*pow（x0,2）/2+x[4]*pow（x1,2）/2-（f[3]-x[3]/6）+（f[4]-x[4]/6）;
printf（"计算结果：
f（4.563）的函数值是：
%lf\nf（4.563）的导数值是：
%lf\n",s4,g4）;}
2.3运行结果：
2.4问题讨论
1.三次样条插值效果比Lagrange插值好，没有Runge现象，光滑性较好。
2.本题的对任意划分的三弯矩插值法可以解决非等距节点的一般性问题。
3.编程过程中由于定义的数组比较多，需要仔细弄清楚各数组所代表的参数，要注意各下标代表的含义，特别是在用追赶法计算的过程中。
3.用Romberg算法求
.
3.1理论依据：
Romberg算法的计算步骤如下：
（1）先求出按梯形公式所得的积分值
（2）把区间2等分，求出两个小梯形面积之和，记为
，即
这样由外推法可得
。
（3）把区间再等分（即
等分），得复化梯形公式
，由
与
外推可得
，
，如此，若已算出
等分的复化梯形公式
，则由Richardson外推法，构造新序列
，m=1,2,…,l,k=1,2,…,l-m+1,
最后求得
。
（4）
或
<
就停止计算，否则回到（3），计算
，一般可用如下算法：
其具体流程如下，并全部存入第一列
通常计算时，用固定l=N来计算，一般l=4或5即能达到要求。
3.2C语言程序代码：
#include
#include
doublef（doublex）//计算f（x）的值
{doublez;
z=pow（3,x）*pow（x,1.4）*（5*x+7）*sin（x*x）;
return（z）;}
main（）
{doublet[20][20],s,e=0.00001,a=1,b=3;
inti,j,l,k;
t[0][1]=（b-a）*（f（b）+f（a））/2;//下为romberg算法
t[1][1]=（b-a）*（f（b）+2*f（（b+a）/2）+f（a））/4;
t[0][2]=（a*t[1][1]-t[0][1]）/（4-1）;j=3;
for（l=2;fabs（t[0][j-1]-t[0][j-2]）>=e;l++）
{for（k=1,s=0;k<=pow（2,l-1）;k++）
s+=f（a+（2*k-1）*（b-a）/pow（2,l））;//判断前后两次所得的T（0）的差是否符合要求，如果符合精度要求则停止循环
t[l][1]=（t[l-1][1]+（b-a）*s/pow（2,l-1））/2;
for（i=l-1,j=2;i>=0;i--,j++）
t[i][j]=（pow（4,j-1）*t[i+1][j-1]-t[i][j-1]）/（pow（4,j-1）-1）;}
if（t[0][1]printf（"t=%0.6f\n",t[0][1]）;
else
printf（"用Romberg算法计算函数所得近似结果为：
\nf（x）=%0.6f\n",t[0][j-1]）;}
3.3运行结果：
3.4MATLAB上机程序
function[T,n]=mromb（f,a,b,eps）
ifnargin<4,eps=1e-6;end
h=b-a;
R（1,1）=（h/2）*（feval（f,a）+feval（f,b））;
n=1;J=0;err=1;
while（err>eps）
J=J+1;h=h/2;S=0;
fori=1:
n
x=a+h*（2*i-1）;
S=S+feval（f,x）;
end
R（J+1,1）=R（J,1）/2+h*S;
fork=1:
J
R（J+1,k+1）=（4^k*R（J+1,k）-R（J,k））/（4^k-1）;
end
err=abs（R（J+1,J+1）-R（J+1,J））;
n=2*n;
end
R;
T=R（J+1,J+1）;
formatlong
f=@（x）（3.^x）*（x.^1.4）*（5*x+7）*sin（x*x）;
[T,n]=mromb（f,1,3,1.e-5）
3.5问题讨论：
1.Romberge算法的优点是:
把积分化为代数运算,而实际上只需求T1（i）,以后用递推可得.算法简单且收敛速度快,一般4或5次即能达到要求。
2.Romberge算法的缺点是:
对函数的光滑性要求较高，计算新分点的值时，这些数值的个数成倍增加。
3.该程序较为复杂，涉及函数定义，有循环，而且循环中又有判断，编写时需要注意该判断条件是处于循环中，当达到要求时跳出循环，终止运算。
4.函数的定义可放在主函数前也可在主程序后面。
本程序采用的后置方式。
4.用定步长四阶Runge-Kutta求解
h=0.0005，打印yi（0.025），yi（0.045），yi（0.085），yi（0.1），（i=1，2，3）
4.1理论依据：
Runge_Kutta采用高阶单步法，这里不是先按Taylor公式展开，而是先写成
处附近的值的线性组合（有待定常数）再按Taylor公式展开，然后确定待定常数，这就是Runge-Kutta法的思想方法。
本题采用四阶古典的Runge-Kutta公式：
4.2C语言程序代码：
#include
voidfun（doublex[4],doubley[4],doubleh）
{y[1]=1*h;
y[2]=x[3]*h;
y[3]=（1000-1000*x[2]-100*x[2]-100*x[3]）*h;//微分方程向量函数}
voidmain（）
{doubleY[5][4],K[5][4],m,z[4],e=0.0005;
doubley[5]={0,0.025,0.045,0.085,0.1};
inti,j,k;
for（i=1;i<=3;i++）
Y[1][i]=0;
for（i=1;i<=4;i++）
for（j=1;j<=3;j++）
K[i][j]=0;
for（k=1;k<=5;k++）
{for（m=y[k-1];m<=y[k];m=m+e）
{for（i=1;i<=3;i++）
z[i]=Y[k][i];
fun（z,K[1],e）;
for（i=1;i<=3;i++）
z[i]=Y[k][i]+e*K[2][i]/2;//依此求K1,K2K3的值
fun（z,K[2],e）;
for（i=1;i<=3;i++）
z[i]=Y[k][i]+e*K[2][i]/2;
fun（z,K[3],e）;
for（i=1;i<=3;i++）
z[i]=Y[k][i]+e*K[3][i];
fun（z,K[4],e）;
for（i=1;i<=3;i++）
Y[k][i]=Y[k][i]+（K[1][i]+2*K[2][i]+2*K[3][i]+K[4][i]）/6;//求Yi[N+1]的值}
if（k!
=5）
for（i=1;i<=3;i++）
Y[k+1][i]=Y[k][i];}
printf（"计算结果:
\n"）;
for（i=1;i<5;i++）
{for（j=1;j<=3;j++）
{printf（"y%d[%4.3f]=%-10.8f,",j,y[i],Y[i][j]）;
if（j==3）
printf（"\n"）;}
printf（"\n"）;}
}
4.3运行结果：
4.4问题讨论：
1.定步长四阶Runge-kutta方法是一种高阶单步法法稳定，精度较高，误差小且程序相对简单，存储量少。
不必求出起始点的函数值，可根据精度的要求修改步长，不会由于起始点的误差造成病态。
2.本程序可以通过修改主程序所调用的函数中的表达式来实现对其它函数的任意初值条件求微分计算。
3.程序中运用了大量的for循环语句，因为该公式中涉及大量的求和，且有不同的函数和对不同的数值求值，编程稍显繁琐。
所以编写过程中一定要注意各循环的次数，以免出错。
5.
用列主元消去法求解Ax=b。
5.1理论依据：
列主元素消元法是在应用Gauss消元法的基础上，凭借长期经验积累提出的，是线性方程组一般解法，目的是为避免在消元计算中使误差的扩大，甚至严重损失了有效数字使数据失真，而在每次初等变换前对矩阵作恰当的调整，以提高Gauss消元法的数字稳定性，进而提高计算所得数据的精确度。
即在每主列中取绝对值最大的元素作主元，再做对应的行交换然后消元求解的办法。
具体做法如下：
将方阵A和向量b写成C=（A，b）。
将C的第1列中第1行的元素与其下面的此列的元素逐一进行比较，找到最大的元素
，将第j行的元素与第1行的元素进行交换，然后通过行变换，将第1列中第2到第n个元素都消成0。
将变换后的矩阵
的第二列中第二行的元素与其下面的此列的元素逐一进行比较，找到最大的元素
，将第k行的元素与第2行的元素进行交换，然后通过行变换，将第2列中第3到第n个元素都消成0。
以此方法将矩阵的左下部分全都消成0后再求解。
最终形式如下：
（A,b）～
5.2C语言程序代码
（1）比较该列的元素的绝对值的大小，将绝对值最大的元素通过行变换使其位于主对角线上；
（2）进行高斯消去法变换，把系数矩阵化成上三角形，然后回代求#include"math.h"
#include"stdio.h"
voidHouseholder（doubleA[9][9]）;
voidexpunction（doubleA[9][9],doubleb[9],doublex[9]）;
voidmain（）
{doubleA[9][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.713847,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}};
doubleb[9]=
{2.1874369,33.992318,-25.173417,0.84671695,1.784317,-86.612343,1.1101230,4.719345,-5.6784392};
doublex[9]={0.0};
inti,j;
Householder（A）;
printf（"\nTheResultsofXare:
\n"）;
expunction（A,b,x）;
for（i=1;i<10;i++）
printf（"X%1d=%f\n",i,x[i-1]）;}
voidHouseholder（doubleA[9][9]）
{doubleq[9],u[9],y[9],s,a,kr;
inti,j,k;
for（i=0;i<7;i++）
{s=0;
for（j=i+1;j<9;j++）
s+=A[j][i]*A[j][i];
s=sqrt（s）;
a=s*s+fabs（A[i+1][i]）*s;
for（j=0;j<9;j++）
{if（j<=i）u[j]=0;
elseif（j==i+1）u[j]=A[j][i]+A[j][i]/fabs（A[j][i]）*s;
elseif（j>i+1）u[j]=A[j][i];}
for（k=0;k<9;k++）
{y[k]=0;
for（j=0;j<9;j++）
y[k]+=A[k][j]*u[j];
y[k]/=a;}
kr=0;
for（k=0;k<9;k++）
kr+=y[k]*u[k];
kr/=2*a;
for（k=0;k<9;k++）q[k]=y[k]-kr*u[k];
for（k=0;k<9;k++）
{for（j=0;j<9;j++）
A[k][j]-=u[k]*q[j]+u[j]*q[k];}
}
}
voidexpunction（doubleA[9][9],doubleb[9],doublex[9]）
{inti,j,k;
doubleB[9][10];
doublez[3];
doublet1=0,t2=0,t3=0;
for（i=0;i<8;i++）
{if（A[i+1][i]>A[i][i]）
{for（j=i,k=0;jz[k]=A[i][j];A[i][j]=A[i+1][j];A[i+1][j]=z[k];
t1=b[i];b[i]=b[i+1];b[i+1]=t1;}
t2=A[i+1][i];
for（j=i;jA[i+1][j]=A[i+1][j]-A[i][j]*t2/A[i][i];
b[i+1]=b[i+1]-b[i]*t2/A[i][i];}
x[8]=b[8]/A[8][8];
for（i=7;i>=0;i--）
{for（j=i+1;j<9;j++）
t3=t3+A[i][j]*x[j];
x[i]=（b[i]-t3）/A[i][i];
t3=0;}
}
5.3运行结果
5.4MATLAB上机程序
unction[x]=mgauss2（A,b,flag）
ifnargin<3,flag=0;end
n=length（b）;
fork=1:
（n-1）
[ap,p]=max（abs（A（k:
n,k）））;
p=p+k-1;
ifp>k
A（[kp],:
）=A（[pk],:
）;
b（[kp],:
）=b（[pk],:
）;
end
m=A（k+1:
n,k）/A（k,k）;
A（k+1:
n,k+1:
n）=A（k+1:
n,k+1:
n）-m*A（k,k+1:
n）;
b（k+
升级会员 