数值分析上机作业参考答案交14789

数值分析上机作业

学院：专业：

学号：姓名

第一题

一、理论依据及算法

1：Household 法

（1） 令A 0=A, a ij (1)

=a ij ，已知A r-1 即A r-1=(a ij (r)

) （2） S r =(

∑

+=n

r i 1

(a ir (r ))2)1/2

（3） αr =S r 2

+|a

(r)

r+1,r |S r

u r =[0,0,a (r)

r+1,r

+Sign(a (r)r+1,r )S r ,a (r)r+2,r ,…,a nr (r)]T

（4） y r =A r-1u r /αr

（5） k r = u r T

y r /2αr （6） q r =y r -k r u r

（7） A r =A r-1-(q r u r T +u r q r T

) r=(1,2,……n-2)

2：超松弛法

其基本思想是在高斯方法已求出x (m)

,x (m-1)

的基础上，组合新的序列，从而

加快收敛速度。其算式是： x i (m)

=(1-ω)x i

(m-1)

+ω(

∑

-=1

1

i j b ij x i (m)

+

∑

+=n

i j 1

x j

(m-1)

+g i )

其中ω是超松弛因子，当ω>1时，可以加快收敛速度。

3：主元消去法：

对矩阵作恰当的调整，选取绝对值尽量大的元素作为主元素。然后把矩阵化为上三角阵，再进行回代，求出方程的解。

二、程序清单：

#include "math.h" main( ) {

int i,j,r,h;

double s, sum,m,m1,n,k,u[9],q[9],y[9],w[9][9], 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.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.0001}}; double y1; r=0;

while(r<7)

{ /*控制循环次数*/

sum=0;

for(h=r+1;h<9;h++) /*求s和n的值*/

sum=sum+A[h][r]*A[h][r];

s=sqrt(sum);

n=s*(s+fabs(A[r+1][r]));

if (A[r+1][r]>=0)

y1=1;

else

y1=-1;

for(i=0;i

u[i]=0;

for(i=r+2;i<9;i++)

u[i]=A[i][r];

u[r+1]=A[r+1][r]+y1*s;

for(j=0;j<9;j++) /*求出y向量*/

{

m=0;

for(i=0;i<9;i++)

m=m+A[i][j]*u[i];

y[j]=m/n;

}

for (i=0;i<9;i++) /*求k和q的值*/

{

m1=0;

for(j=0;j<9;j++)

m1=m1+u[j]*y[j];

k=m1/(2*n);

q[i]=y[i]-k*u[i];

}

for(i=0;i<9;i++) /*求结果*/

for(j=0;j<9;j++)

{

w[i][j]=q[i]*u[j]+u[i]*q[j];

A[i][j]=A[i][j]-w[i][j];

}

r++;

}

printf("转化后的矩阵为：B[9][9]=:\n"); /*打印转化后的矩阵*/ for(i=0;i<=8;i++)

{

for(j=0;j<=8;j++)

printf("%4.4f ",A[i][j]);

printf("\n");

}

}

用超松弛法解方程组

#include"math.h"

main( )

{

int i,j,k,r,m;

double w=1.4,g1[9],B1[9][9],b[9][9],x[9]={0},

g[9]={2.1874369,33.992318,-25.173417,0.84671695,1.784317,

-86.612343,1.1101230,4.719345,-5.6784392},

B[9][9]={{12.384120,-4.893077},{-4.893077,25.398416,6.494097},

{0,6.494097,20.611499,8.243925},{0,0,8.243925,23.422838,-13.880071},

{0,0,0,-13.880071,29.698278,4.534502},{0,0,0,0,4.534502,16.00617,4.881435},

{0,0,0,0,0,4.881435,26.013315,-4.503635},{0,0,0,0,0,0,-4.503635,21.254061,4.504498}, {0,0,0,0,0,0,0,4.504498,14.534122}},

e[9][9]={{1},{0,1},{0,0,1},{0,0,0,1},{0,0,0,0,1},{0,0,0,0,0,1},{0,0,0,0,0,0,1},

{0,0,0,0,0,0,0,1},{0,0,0,0,0,0,0,0,1}};

i=0;

j=0;

while(i<9) /*求出矩阵b和g1的值*/

{

while(j<9)

{

B1[i][j]=B[i][j]/B[i][i];

b[i][j]=e[i][j]-B1[i][j];

j++;

}

g1[i]=g[i]/B[i][i];

i++;

}

k=0;

while(k<9) /*执行本算法*/

{

x[0]=(1-w)*x[0]+w*(b[0][1]*x[1]+g1[0]);

for(r=1;r<8;r++)

x[r]=(1-w)*x[r]+w*(b[r][r-1]*x[r-1]+b[r][r+1]*x[r+1]+g1[r]);

x[8]=(1-w)*x[8]+w*(b[8][7]*x[7]+g1[8]);

k++;

}

for(m=0;m<9;m++) /*输出结果*/

printf("%f\n",x[m]);

}

用消去法解方程组

#include"math.h"

main( )

{

int i,j,k,n;

double x[9],m[9],l[9],h[9],

b[9]={2.1874369,33.992318,-25.173417,0.84671695,1.784317,-86.612343,

1.1101230,4.719345,-5.6784392},

B[9][9]={{12.3841,-4.8931,},{-4.8931,25.3984,6.4941,},{0,6.4941,20.6116,8.2440,},

{0,0,8.2440,23.4231,-13.8802,},{0,0,0,-13.8802,29.6980,4.5344,},{0,0,0,0,4.5344,16.0061, 4.8814,},{0,0,0,0,0,4.8814,26.0133,-4.5036,},{0,0,0,0,0,0,-4.5036,21.2541,4.5045},{0,0,0,0,0, 0,0,4.5045,14.5341}};

i=0;

while(i<8) /*迭代次数*/

{

l[i]=B[i+1][i]/B[i][i];

for(j=0;j<8;j++) /*求矩阵B和b的值*/

{

h[j]=B[i][j]*l[i];

B[i+1][j]=B[i+1][j]-h[j];

}

b[i+1]=b[i+1]-b[i]*l[i];

i++;

}

k=7;

while(k>0) /*求出x的值*/

{

x[8]=b[8]/B[8][8];

m[k]=b[k]-B[k][k+1]*x[k+1];

x[k]=m[k]/B[k][k];

m[0]=b[0]-B[0][1]*x[1];

x[0]=m[0]/B[0][0];

k--;

}

for(n=0;n<9;n++)

printf("%f\n",x[n]);

}

三、运行结果：

第一题（1）结果为：B[9][9]=

12.3841 –4.8931 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000

-4.8931 25.3984 6.4941 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000

0.0000 6.4941 20.6116 8.2440 0.0000 0.0000 0.0000 0.0000 0.0000

0.0000 0.0000 8.2440 23.4231 -13.8802 0.0000 0.0000 0.0000 0.0000

0.0000 0.0000 0.0000 -13.8802 29.6980 4.5344 0.0000 0 .000 0.0000

0.0000 0.0000 0.0000 0.0000 4.5344 16.0061 4.8814 0.0000 0.0000

0.0000 0.0000 0.0000 0.0000 0.0000 4.8814 26.0133 -4.5036 0.0000

0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 -4.5036 21.2541 4.5045

0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 4.5045 14.5341

第一题（2） 结果为：

x[0]=1.073409 x[1]=2.272579 x[2]=-2.856599 x[3]=2.292510 x[4]=2.112159 x[5]=-6.422562 x[6]=1.357798 x[7]=0.634258 x[8]=-0.587041 第一题（3） 结果为：

x[0]=1.075802 x[1]=2.275736 x[2]=–2.855472 x[3]=2.293009 x[4]=2.112517 x[5]=–6.423293 x[6]= 1.356278 x[7]=0.625385 x[8]=–0.547122

四、问题讨论

如果所解方程系数矩阵B 是三对角阵，SOR 方法的矩阵形式为： X (m)=(E-ωL)-1((1-ω)E+ωR)x (m-1)+(E-ωL)-1

ωg

若记 L ω=(E-ωL)-1

((1-ω)E+ωR)，SOR 收敛的充要条件是S(L ω)<1.且若A 为对称正定阵，则当松弛因子ω满足0<ω<2时，SOR 方法收敛。用SOR 方法解时，收敛速度很快。用消去法解时，可以用特殊的解法，效果较好。

第三题

1． 题目： 用三次样条插值求)563

.4(f 及)563.4(f '的近似值。 2． 算法：

∑∑∑-==-=-+Ω=-Ω=-Ω

=

101

12

1

3310

1

3

)2()()(

)(j j j j j j

j

j x x j x c h

x x c x s

∑+=++-++

-=

Ω10

1!)2

1

()1()(k j k j k j k k j k x C x 其中j c 必须的方程 ：

???

??

?

??

?????

?

??

?

?

?

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

??---01551.1481551.1318338.1347667.12675460.11750557.106566268.93177664.85916738.6158839.402241411411411411411

41141141141141421098765432101c c c c c c c c c c c c

q[0]=0 , u[0]=0 ,

1,,2,1]),[][][(][][-???=?+-=n i i q i a i b i c i q

n i i q i a i b i u i a i d i u ,,2,1]),1[][][(])[][][(}[???=-?+?-=

x[9]=u[9]

1,,2,1],[]1[][][???--=++?=n n i i u i x i q i x []

3

33333)2()1(46)1(4)2(6

1)(+++++-+--++-+=

Ωx x x x x x s '(x)= ∑∑-=-=--Ω--+Ω10

1

210

12)21

()21(j j j j j x c j x c

2

)3(2

4c c s -=

' 3．程序清单：

#include "math.h" #include "stdio.h" #define M 12

void zhuigan(double coeff[M][M],double b[M],double c[M],int n) {

int i,j,m=n+3; double s1,s2,q[M],u[M+1]; s1=coeff[0][0];s2=b[0]; i=1;

do {if(fabs(s1)<=1e-6){printf("\nThe diagral element is zero!");return;} else {q[i]=-coeff[i-1][i]/s1; u[i]=s2/s1;} i++;

s1=coeff[i-1][i-1]+coeff[i-1][i-2]*q[i-1]; s2=b[i-1]-coeff[i-1][i-2]*u[i-1];} while(i

if(fabs(s1)<=1e-6){printf("\nThe yangtiao element is zero!"); return;} else u[m]=s2/s1; c[m-1]=u[m];

for(i=m-2;i>=0;i--) c[i]=q[i+1]*c[i+1]+u[i+1];} double Bas_3(double x)

{double t,a,y; t=fabs(x); a=x*x; if(t>=2.0) y=0;

else if(t>1.0&&t<2.0) y=-a*t/6+a-2*t+4.0/3; else y=a*t/2-a+2.0/3; return y;}

double Bas_2(double x) {double y,t; t=fabs(x); if(t>=3.0/2) y=0;

else if(t>1.0/2&&t<3.0/2) y=t*t/2-3*t/2+9.0/8; else y=-t*t+3.0/4; return y;}

double DBas_3(double x)

{return (Bas_2(x+1.0/2)-Bas_2(x-1.0/2));} main()

{int i,j,n,m; double h,t,x,dx,s,ds,x0, coeff[M][M],b[M],c[M],

a[10]={0,0.69314718,1.0986123,

1.3862944,1.6094378,1.7917595,1.9459101,

2.079445,2.1972246,2.3025851},da[2]={ 1,0.1};

printf("\n\n\n ");

printf("\n输入初始等分区间个数:n="); scanf("%d",&n); /*n+3为插值条件数*/ printf("\n输入步长:h="); scanf("%lf",&h);

printf("\n输入函数自变量起始值:x0=");scanf("%lf",&x0);

coeff[0][0]=2;coeff[0][1]=4; coeff[n+2][n+2]=2;coeff[n+2][n+1]=4;

for(i=1;i

b[0]=da[0]*2*h; b[n+2]=da[1]*2*h;

for(i=1;i

b[0]=b[1]-b[0];b[n+2]=b[n+1]+b[n+2];

zhuigan(coeff,b,c,n);

printf("\n输入f(x)插值点:x="); scanf("%lf",&x);

for(j=0,s=0;j

printf("\n输入f'(x)插值点:dx="); scanf("%lf",&dx);

for(j=0,ds=0;j

ds/=h;

printf("\n结果为:f(%.4f)=%.8f\nf'(%.4f)=%.8f\n",x,s,dx,ds);}

运行结果：

输入初始等分区间个数:n=9

输入步长:h=1

输入函数自变量起始值:x0=1

输入f(x)插值点:x=4.563

输入f'(x)插值点:dx=4.563

结果为: f(4.5630)=1.51793242 f'(4.5630)=0.21928608

5．问题讨论

样条插值效果比Lagrange插值好,没有Runge现象.

第四题

一、理论依据及算法

设函数在有限区间[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迭代过程为

x k+1=φ(x k)=x-f(x k)/f’(x k),k=1,2,3…

令

)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

二、程序清单

#include

double f(double x) /*定义函数f(x)*/ {

return (x*x*x*x*(x*x*x-28)+14); }

double f1(double x) /*定义一阶导函数*/ {

return (x*x*x*(7*x*x*x-28*4)); }

double f2(double x) /*定义二阶导函数*/ {

return (x*x*(42*x*x*x-28*4*3)); }

main() {

double x0,x1;

if (f(1.9)*f2(1.9)>0) /*判断条件*/ x0=1.9;

else x0=0.1; do {

x1=x0;

x0=x0-f(x0)/f1(x0); }

while (fabs(x1-x0)>=1.0e-5); /*设定停机条件*/ printf("x=%f",x0); }

三、 运行结果

x=0.845497

第五题

一 程序要求 用Romberg 算法求

?=+3

1

24.1)00001.0(,sin )75(3

εdx x x x x

二 程序算法

引入记号T m (k)称为Romberge 序列,当|T i (0)-T i+1(0)

|< 0.00001就停止计算:

T 1(0)

=(b-a)(f(a)-f(b))/2

T 1(l)

=1/2[T 1

(l-1)

+b-a/2

l-1

∑

=2

1

i f(a+(2i-1)b-a/2)]

T m+1(k-1)

=4m

T m (k)-T m (k-1)

/4m

-1

三 程序

#define N 12

#include "math.h" #include "stdio.h"

#define Epxilun 0.00001 /*允许误差*/ float f(float x) {

float sum=0.0;

sum=pow(3,x)*pow(x,1.4)*sin(x*x)*(5*x+7); return(sum); }

main() {

float T[N][N];

float a=1,b=3,sum=0.0; int j,i;

T[0][0]=(b-a)*(f(a)+f(b))/2; j=1; do {

sum=0.0;

for(i=1;i<=pow(2,j-1);i++)

sum=sum+f(a+(2*i-1)*(b-a)/pow(2,j)); sum=sum*(b-a)/pow(2,j-1); T[j][0]=(T[j-1][0]+sum)/2; for(i=1;i<=j;i++)

T[j-i][i]=(pow(4,i)*T[j-i+1][i-1]-T[j-i][i-1])/(pow(4,i)-1); j++; }

while(fabs(T[0][j-1]-T[0][j-2])>=Epxilun);

printf("\n THE romberg result is: T0=%f T1=%f ",T[0][j-1],T[0][j-2]); getch(); }

四 运行结果：

T0=440.536011, T1=440.536011

∫313xx1.4(5x+7)sinx2dx=440.536011

五 问题讨论

Romberge 算法的优点是:

1. 把积分化为代数运算,而实际上只需求T1(i),以后用递推可得.

2. 算法简单且收敛速度快,一般4或5次即能达到要求.

3. 节省存储量.

Romberge 算法的缺点是:

1. 对函数的光滑性要求较高.

第六题

1． 用定步长四Runge-Kutta 求解常微分方程组。

1/1=dt dy

32/y dt dy =

32310010001000/y y dt dy --= ()001=y ()002=y

()003=y

()()()()()3,2,1,1.0,085.0,045.0,025.0,0005.0==i y y y y h i i i 打印

2． 算法：

][j y i 表示3,2,1),(=i hj y i

????

?

??????

++=++=++==++++=+])[][,(][])

[21][,21(][])

[21][,21(][][])[][2][2][(6

1][]1[34231214321i hk n y h t f i k i hk n y h t f i k i hk n y h t f i k f i k i k i k i k i k n y n y i in i in i in in i i

3．程序清单：

#include "stdio.h"

void F(double Y[3],double K[3])

{

K[0]=1;

K[1]=Y[2];

K[2]=1000-1000*Y[1]-100*Y[2];

}

main()

{

double Y[4][3],h=0.0005,m,K1[3],K2[3],K3[3],K4[3],z[3];

double y[4]={0.025,0.045,0.085,0.1};

int i,j,k;

for(i=0;i<3;i++)

Y[0][i]=0;

m=0;

k=0;

while(k<5)

{

for(;m

{

for(i=0;i<3;i++)

z[i]=Y[k][i];

F(z,K1);

for(i=0;i<3;i++)

z[i]=Y[k][i]+0.5*h*K1[i];

F(z,K2);

for(i=0;i<3;i++)

z[i]=Y[k][i]+0.5*h*K2[i];

F(z,K3);

for(i=0;i<3;i++)

z[i]=Y[k][i]+h*K3[i];

F(z,K4);

for(i=0;i<3;i++)

Y[k][i]+=h/6*(K1[i]+2*K2[i]+2*K3[i]+K4[i]);

}

if(k!=4)

for(j=0;j<3;j++)

Y[k+1][j]=Y[k][j];

k++;

}

printf("结果为:\n");

for(i=0;i<4;i++)

{

for(j=0;j<3;j++)

printf("y%d(%4.3f)=%f\n",j+1,y[i],Y[i][j]);

printf("\n");

}

}

计算结果如下：

y1(0.025)=0.025000

y2(0.025)=0.151599

y2(0.025)=8.335426

y1(0.045)=0.045000

y2(0.045)=0.312860

y3(0.045)=7.536280

y1(0.085)=0.085000

y2(0.085)=0.560581

y3(0.085)=4.946353

y1(0.100)= 0.10000

y2(0.100)=0.628881

y3(0.100)=4.180993

5．问题讨论

Runge_Kutta方的优点:

1.精度高,不必用别的方法求开始几点的函数值.

2.可根据f’(t,y)变化的情况与需要的精度自动修改步长.

3.程序简单,存储量少.

4.方法稳定.

Runge_Kutta方的缺点:

每步要计算函数值f(t,y)四次,在f(t,y)，较复杂时,工作量大, 且每一步缺乏可靠的检查.