《数值分析》上机实验报告.doc

数值分析上机实验报告

《数值分析》上机实验报告
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.2 C语言程序原代码：
#include<stdio.h>
#include<math.h>
main()
{double x2,f,f1;
double x1=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.4 MATLAB上机程序
function y=Newton(f,df,x0,eps,M)
d=0;
for k=1:M
if feval(df,x0)==0
d=2;break
else
x1=x0-feval(f,x0)/feval(df,x0);
end
e=abs(x1-x0);
x0=x1;
if e<=eps&&abs(feval(f,x1))<=eps
d=1;break
end
end
if d==1
y=x1;
elseif d==0
y='迭代M次失败';
else
y= '奇异'
end
function y=df(x)
y=7*x^6-28*4*x^3;
End
function y=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.2 C语言程序代码：
#include<stdio.h>
#include<math.h>

void main()
{int i,j,m,n,k,p;
double q10,p10,s4,g4,x0,x1,g0=1,g9=0.1;;
double s[10][10];
double a[10],b[10],c[10],d[10],e[10],x[10],h[9],u[9],r[9];
double f[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.2 C语言程序代码：
#include<math.h>
#include<stdio.h>
double f(double x) //计算f(x)的值
{double z;
z=pow(3,x)*pow(x,1.4)*(5*x+7)*sin(x*x);
return(z);}
main()
{ double t[20][20],s,e=0.00001,a=1,b=3;
int i,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]<e)
printf("t=%0.6f\n",t[0][1]);
else
printf("用Romberg算法计算函数所得近似结果为：\nf(x)=%0.6f\n",t[0][j-1]);}

3.3 运行结果：

3.4 MATLAB上机程序
function [T,n]=mromb(f,a,b,eps)
if nargin<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;
for i=1:n
x=a+h*(2*i-1);
S=S+feval(f,x);
end
R(J+1,1)=R(J,1)/2+h*S;
for k=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);

format long
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.2 C语言程序代码：
#include<stdio.h>
void fun(double x[4],double y[4],double h)
{y[1]=1*h;
y[2]=x[3]*h;
y[3]=(1000-1000*x[2]-100*x[2]-100*x[3])*h; //微分方程向量函数}
void main()
{ double Y[5][4],K[5][4],m,z[4],e=0.0005;
double y[5]={0,0.025,0.045,0.085,0.1};
int i,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.2 C语言程序代码
（1）比较该列的元素的绝对值的大小，将绝对值最大的元素通过行变换使其位于主对角线上；
（2）进行高斯消去法变换，把系数矩阵化成上三角形，然后回代求#include "math.h"
#include "stdio.h"
void Householder(double A[9][9]);
void expunction(double A[9][9],double b[9],double x[9]);
void main()
{double A[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}};
double b[9]=
{2.1874369,33.992318,-25.173417,0.84671695,1.784317,-86.612343,1.1101230,4.719345,-5.6784392};
double x[9]={0.0};
int i,j;
Householder(A);
printf("\n The Results of X are:\n");
expunction(A,b,x);
for(i=1;i<10;i++)
printf("X%1d=%f\n",i,x[i-1]);}
void Householder(double A[9][9])
{double q[9],u[9],y[9],s,a,kr;
int i,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;
else if(j==i+1) u[j]=A[j][i]+A[j][i]/fabs(A[j][i])*s;
else if(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];}
}
}
void expunction(double A[9][9],double b[9],double x[9])
{int i,j,k;
double B[9][10];
double z[3];
double t1=0,t2=0,t3=0;
for(i=0;i<8;i++)
{if(A[i+1][i]>A[i][i])
{for(j=i,k=0;j<i+3;j++,k++)
z[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;j<i+3;j++)
A[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.4 MATLAB上机程序
unction [x]=mgauss2(A,b,flag)
if nargin<3,flag=0;end
n=length(b);
for k=1:(n-1)
[ap,p]=max(abs(A(k:n,k)));
p=p+k-1;
if p>k
A([k p],:)=A([p k],:);
b([k p],:)=b([p k],:);
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+1:n)=b(k+1:n)-m*b(k);
A(k+1:n,k)=zeros(n-k,1);
if flag~=0,Ab=[A,b],end
end
x=zeros(n,1);
x(n)=b(n)/A(n,n);
for k=n-1:-1:1
x(k)=(b(k)-A(k,k+1:n)*x(k+1:n))/A(k,k);
end

format long
A=[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];
b=[2.1874369,33.992318,-25.173417,0.84671695,1.784317,-86.612343,1.1101230,4.719345,-5.6784392]';
x=mgauss2(A,b);x=x'

5.5问题讨论：
1.输入矩阵可用循环比较指令确认输入的是否为对称矩阵。
2.循环体中的累加值注意初始化零。
3.注意公式中的下标从一开始，数组中的下标从零开始。
4.在编程中数组的角标从0开始与数学中的起始脚表不同，编程时必须注意。
5.在给数组元素通过表达式赋值时要注意原始赋值的覆盖问题。其中专业理论知识内容包括：保安理论知识、消防业务知识、职业道德、法律常识、保安礼仪、救护知识。作技能训练内容包括：岗位操作指引、勤务技能、消防技能、军事技能。
二．培训的及要求培训目的
安全生产目标责任书
为了进一步落实安全生产责任制，做到“责、权、利”相结合，根据我公司2015年度安全生产目标的内容，现与财务部签订如下安全生产目标：
一、目标值：
1、全年人身死亡事故为零，重伤事故为零，轻伤人数为零。
2、现金安全保管，不发生盗窃事故。
3、每月足额提取安全生产费用，保障安全生产投入资金的到位。
4、安全培训合格率为100%。
二、本单位安全工作上必须做到以下内容：
1、对本单位的安全生产负直接领导责任，必须模范遵守公司的各项安全管理制度，不发布与公司安全管理制度相抵触的指令，严格履行本人的安全职责，确保安全责任制在本单位全面落实，并全力支持安全工作。
2、保证公司各项安全管理制度和管理办法在本单位内全面实施，并自觉接受公司安全部门的监督和管理。
3、在确保安全的前提下组织生产，始终把安全工作放在首位，当“安全与交货期、质量”发生矛盾时，坚持安全第一的原则。
4、参加生产碰头会时，首先汇报本单位的安全生产情况和安全问题落实情况；在安排本单位生产任务时，必须安排安全工作内容，并写入记录。
5、在公司及政府的安全检查中杜绝各类违章现象。
6、组织本部门积极参加安全检查，做到有检查、有整改，记录全。
7、以身作则，不违章指挥、不违章操作。对发现的各类违章现象负有查禁的责任，同时要予以查处。
8、虚心接受员工提出的问题，杜绝不接受或盲目指挥；
9、发生事故，应立即报告主管领导，按照“四不放过”的原则召开事故分析会，提出整改措施和对责任者的处理意见，并填写事故登记表，严禁隐瞒不报或降低对责任者的处罚标准。
10、必须按规定对单位员工进行培训和新员工上岗教育；
11、严格执行公司安全生产十六项禁令，保证本单位所有人员不违章作业。
三、 安全奖惩：
1、对于全年实现安全目标的按照公司生产现场管理规定和工作说明书进行考核奖励；对于未实现安全目标的按照公司规定进行处罚。
2、每月接受主管领导指派人员对安全生产责任状的落

21