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

数值分析上机实验报告

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

1.用Newton 法求方程 X 7-X 4+14=0

在（0.1，1.9）中的近似根（初始近似值取为区间端点，迭代6次或误差小于0.00001）。 1.1 理论依据：

设函数在有限区间[a ，b]上二阶导数存在，且满足条件

{}α?上的惟一解在区间平方收敛于方程所生的迭代序列

迭代过程由则对任意初始近似值达到的一个中使是其中上不变号在区间],[0)(3,2,1,0,)

(')

()(],,[x |))(),((|,|,)(||

)(|.

4;

0)(.3],[)(.20

)()(.110......b a x f x k x f x f x x x Newton b a b f a f m ir b a c x f a

b c f x f b a x f b f x f k k k k k k ==-

==∈≤-≠>+

令

)9.1()9.1(0)8(4233642)(0)16(71127)(0)9.1(,0)1.0(,1428)(3

2

2

5

333647>?''<-=-=''<-=-='<>+-=f f x x x x x f x x x x x f f f x x x f

故以1.9为起点

??

??

?

='-

=+9.1)()(01x x f x f x x k k k k 如此一次一次的迭代，逼近x 的真实根。当前后两个的差<=ε时，就认为求出了近似的根。本程序用Newton 法求代数方程(最高次数不大于10)在（a,b ）区间的根。

1.2 C语言程序原代码：

#include

#include

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.已知函数值如下表：

试用三次样条插值求f(4.563)及f ’(4.563)的近似值。 2.1 理论依据

33

22111

11

111

1

11

()()()()()()()

6666j j j j j j j j

j j j j j j j j x x x x h x x h x x S x M M y M y M h h h h ---------------=++-+-这里11j j j h x x --=- ，所以只要求出j M ，就能得出插值函数S （x ）。

求j M 的方法为：001111221

121

22

212N N N N M d M d M d μ

λμλμλ--??????

??????????????????=?

???????????

???????????????????

?????

这里1000

00

11111

111

116()6()(1,2,,1)61[()]

1j j j j j

j j j j N N N N N N j j

j j j j j j j y y d y h h

y y y y d j N h h h h d y y y h h h h h h h h μλμ+----------?

'=-???--=-=-?+?

?

?'=--??

?==-=?++?

最终归结为求解一个三对角阵的解。

用追赶法解三对角阵的方法如下：

1

1112221222111111111n n n n n n n n n b c a b c l A LU l b c a l a b γβγββγβ-----?????????????

???????????===?

??????

?????????????????

1,,,n L d LUx d L d Ux δδδδδδ??

=???===???=?????

即若记则由得 112111n n n d l l d δδ??????

????????????=????????????

?????? ， 1

1

111

1n n n n n x x βγδβγβδ--??????????????????=??????

???????

?????

综上可得求解方程Ax=d 的算法：

111

1111111111,,,,1,2,3,,1,,1,,2,1

i i i i i i i

i i i i

n i i i n i n i b d l b l c d l i n c x

x x i n αβδββδδδδββ+++++++++?

====-???

=-=-??-?===-??

2.2 C 语言程序代码：

#include

#include

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 算法求)00001.0(sin )75(323

14.1=+?ε允许误差dx x x x x . 3.1 理论依据：

Romberg 算法的计算步骤如下：

（1）先求出按梯形公式所得的积分值

(0)1[()()]2

b a

T f a f b -=

+ （2）把区间2等分，求出两个小梯形面积之和，记为(1)1T ，即

(1)1[()()2()]42

b a b a

T f a f b f -+=

++ 这样由外推法可得(0)2T ,(1)2(0)

(1)(0)11(0)

11221

()421411()2

T T T T T --==

--。 （3）把区间再等分（即22等分），得复化梯形公式(2)1T ，由(1)1T 与(2)1T 外推可得(2)(1)(1)112

441T T T

-=-，2(1)(0)(0)2232

441

T T T -=-，如此，若已算出2k 等分的复化梯形公式()1k T ，则由Richardson 外推法，构造新序列 ()(1)(1)

1

441

m k k k m m

m m T T T

--+-=

-， m=1,2,…,l, k=1,2,…,l-m+1, 最后求得(0)1l T +。

（4）(0)(0)1l l T T +≈或(0)(0)

1

||l l T T +-<ε就停止计算，否则回到（3），计算(1)1l T +，一般可用如下算法：

1(0)12()(1)

1111()(1)(1)

1[()()]21{[(21)]}2224,1,2,,,1,2,,141l l l l l

i m k k k m m m m

b a

T f a f b b a b a T T f a i T T T m l k l m ---=--+?-=+??

?--?=++-??

?-?===-+-??

∑

其具体流程如下，并全部存入第一列

(0)

1(1)T 2(3)T (0)3(6)T (1)1(2)T 2(5)T (2)1(4)T 通常计算时，用固定l=N 来计算，一般l=4或5即能达到要求。

3.2 C 语言程序代码：

#include #include

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]

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 求解

??????????

?===--===0

)0(0)0(0)0(10010001000//1/321

3

2332

1y y y y y dt dy y

dt dy dt dy h =0.0005，打印y i (0.025) ， y i (0.045) ， y i (0.085) ， y i (0.1) ，(i =1，2，3) 4.1 理论依据：

Runge_Kutta 采用高阶单步法，这里不是先按Taylor 公式展开，而是先写成n t 处附近的值的线性组合（有待定常数）再按Taylor 公式展开，然后确定待定常数，这就是Runge-Kutta 法的思想方法。

本题采用四阶古典的Runge-Kutta 公式：

)

,()3/,3/2()3/,3/()

,(8/]33[321421312143211hK hK hK Y h x hF K hK hK Y h x hF K hK Y h x hF K Y x hF K K K K K Y Y n n n n n n n n n n +-++=+++=++==++++=+

4.2 C 语言程序代码：

#include

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.

??

?

?????

?????

???????????????=40.00001 4.446782 2.213474- 0.238417 1.784317 0.037585- 1.010103- 3.123124 2.031743- 4.446782 30.719334 3.123789 1.103456- 2.121314 0.71828- 0.336993 1.112348 3.067813 2.213474- 3.123789 14.7138465 0.103458- 3.111223- 2.101023 1.056781- 0.784165- 1.7423820.238417 1.103456- 0.103458- 9.789365 0.431637 3.741856- 1.836742 1.563849 0.718719 1.784317

2.121314

3.111223- 0.431637 19.897918

4.101011 2.031454 2.189736 0.113584-0.037585- 0.71828- 2.101023 3.741856- 4.101011 27.108437 3.123848 1.012345- 1.112336 1.010103- 0.336993 1.056781- 1.836742 2.031454 3.123848 1

5.567914 3.125432- 1.061074- 3.123124 1.112348 0.784165- 1.563849 2.189736 1.012345- 3.125432- 19.141823 2.115237 2.031743- 3.067813 1.742382 0.718719 0.113584- 1.112336 1.061074- 2.115237 12.38412

A T

b )5.6784392- 4.719345 1.1101230 86.612343- 1.784317 0.84671695 25.173417- 33.992318 2.1874369(= 用列主元消去法求解Ax=b 。

5.1 理论依据：

列主元素消元法是在应用Gauss 消元法的基础上，凭借长期经验积累提出的，是线性方程组一般解法，目的是为避免在消元计算中使误差的扩大，甚至严重损失了有效数字使数据失真，而在每次初等变换前对矩阵作恰当的调整，以提高Gauss 消元法的数字稳定性，进而提高计算所得数据的精确度。即在每主列中取绝对值最大的元素作主元，再做对应的行交换然后消元求解的办法。具体做法如下：

将方阵A 和向量b 写成C=（A ，b ）。将C 的第1列中第1行的元素与其下面的此列的元素逐一进行比较，找到最大的元素1j c ，将第j 行的元素与第1行的元素进行交换，然后通过行变换，将第1列中第2到第n 个元素都消成0。将变换后的矩阵(1)C 的第二列中第二行的元素与其下面的此列的元素逐一进行比较，找到最大的元素(1)2k c ，将第k 行的元素与第2行的元素进行交换，然后通过行变换，将第2列中第3到第n 个元素都消成0。以此方法将矩阵的左下部分全都消

成0后再求解。最终形式如下：

(A,b）～

()()

1 111

()

00

n n

n

n

n

nn

g a a

g

a

?? ?

*

? ?

?

??

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.4 46782},

{-2.031743,3.123124,-1.010103,-0.037585,1.784317,0.238417,-2.213474,4.446782,4 0.00001}};

double b[9]=

{2.1874369,33.992318,-25.173417,0.84671695,1.784317,-86.612343,1.1101230,4.71 9345,-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

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

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);