东南大学数值分析上机作业

习题117．〔上机题〕舍入误差与有效数设，其精确值为。〔1〕编制按从大到小的顺序，计算的通用程序。〔2〕编制按从小到大的顺序，计算的通用程序。〔3〕按两种顺序分别计算，，，并指出有效位数。〔编制程序时用单精度〕〔4〕通过本上机题，你明白了什么？答：1、从大到小顺序：#include

<stdio.h>

int

main()

{

float

S=0,T,N,i=2;

printf("请输入N：");

scanf("%f",&N);

while

(i<=N)

{

T=1/(i*i-1);

i=++i;

S+=T;

}

printf("S=%f

\n",S);

system("pause");

return

0;

}

2、从小到大顺序：#include

<stdio.h>

int

main()

{

float

S=0,T,N;

printf("请输入N：");

scanf("%f",&N);

while

(N>=2)

{

T=1/(N*N-1);

N=--N;

S+=T;

}

printf("S=%f

\n

",S);

system("pause");

return

0;

}

3、计算结果：N101010从大到小S0.7400490.7498520.749852从小到大S0.7400500.7499000.749999精确值0.7400490.74990.749999有效位数从大到小S643从小到大S5464、感想：算法的不同导致计算结果的误差不同，一定要好好设计算法，尽量使算法精确。通过上面的例子，可以看出按不同的计算顺序所得到的结果是不同的，按从大到小的顺序计算的值与精确值有较大的误差，而按从小到大的顺序计算的值与精确值吻合。计算机在进行数值计算时会出现“大数吃小数〞的现象，导致计算结果的精度有所降低，使用计算机进行同号数的加法时，采用绝对值较小者先加的算法，得到的结果的相对误差较小。5、附图：图SEQ图\*ARABIC1从大到小算法及结果图SEQ图\*ARABIC2从小到大算法及结果习题220．〔上机题〕Newton迭代法〔1〕给定初值及容许误差，编制Newton法解方程根的通用程序。〔2〕给定方程，易知其有三个根，，。1．由Newton方法的局部收敛性可知存在，当时，Newton迭代序列收敛于根。试确定尽可能大的。2．试取假设干初始值，观察当，，，，时Newton序列是否收敛以及收敛于哪一个根。〔3〕通过本上机题，你明白了什么？答：1、通用程序：#include<stdio.h>

#include<math.h>

float

f(float

x)

{

float

f;

f=x*x*x/3-x;

return(f);

}

float

df(float

x)

{

float

df;

df=x*x-1;

return

(df);

}

int

main()

{

float

x0,x1,d,eps=0.0000005;

int

k=0;

printf("请输入初值：");

scanf("%f",&x0);

do

{

d=-f(x0)/df(x0);

x1=x0+d;

k++;

x0=x1;

}

while(fabs(d)>eps);

printf("迭代%d次后,根=%f\n",k,x0);

system("pause");

return

0;

}

2、确定尽可能大的：〔由于函数是奇函数，故仅需在正半轴求〕#include<stdio.h>

#include<math.h>

float

f(float

x)

{

float

f;

f=x*x*x/3-x;

return(f);

}

float

df(float

x)

{

float

df;

df=x*x-1;

return

(df);

}

int

main()

{

float

x0=0,x1,d;

float

eps=0.000005,delta=0;

while(x0==0)

{

delta=delta+0.000001;

x0=x0+delta;

do

{

d=-f(x0)/df(x0);

x1=x0+d;

x0=x1;

}

while(fabs(d)>eps);

}

delta=delta-0.000001;

printf("尽可能大的值=%f\n",delta);

system("pause");

return

0;

}

运算结果为：0.7745963、测试是否在给定区间收敛：在五个区间内分别取初值-100000，-0.8，0.5，0.9，110000测试是否收敛及收敛到哪个根，结果如下表：初值-100000-0.80.50.9110000是否收敛是是是是是收敛于-3-30334、感想：选取的迭代初值不同，则迭代序列将有可能收敛于不同的根，牛顿法对迭代初值的要求较高。有些区间收敛范围相当大，有些却非常小。另外，我发现容许误差的大小对迭代结果的精度有很大的影响，如果想得到比拟精确的结果，要将容许误差设置小一点。5、附图：图1利用牛顿法求方程的近似解图2确定收敛范围习题339．〔上机题〕列主元Gauss消去法对于某电路的分析，归结为求解线性方程组RI=V，其中R为9阶矩阵，参见课本127页VT=〔1〕编制解n阶线性方程组Ax=b的列主元三角分解法的通用程序；〔2〕用所编制的程序解线性方程组RI=V，并打印出解向量，保存五位有效数；〔3〕本编程之中，你提高了哪些编程能力？答:1、通用程序#

include<stdio.h>

#

include<math.h>

#

define

delta

1e-6

#define

N

100

int

main()

{

int

i,j,t,r,n,u,c=0;

float

p,L,max,s;

float

X[N];

float

a[N][N+1];

printf("请输入系数矩阵的阶数：\n");

scanf("%d",&n);

printf("输入的增广矩阵系数,中间用空格隔开：\n");

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

for(j=0;j<n+1;j++)

scanf("%f",&a[i][j]);

printf("你输入的增广矩阵为：\n");

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

{

for(j=0;j<n+1;j++)

printf("%f\t",a[i][j]);

printf("\n");

}

for(j=0;j<n-1;j++)

{

{

max=fabs(a[j][j]);

r=j;

}

for(i=j+1;i<n;i++)

if(fabs(a[i][j])>max)

{

max=a[i][j];

r=i;

}

if(fabs(a[i][j])<delta)

printf("矩阵奇异");

for(t=j;t<n+1;t++)

{

p=a[j][t];

a[j][t]=a[r][t];

a[r][t]=p;

}

for(i=j+1;i<N;i++)

{

L=a[i][j]/a[j][j];

for(t=j;t<n+1;t++)

a[i][t]=a[i][t]-L*a[j][t];

}

}

printf("输出原方程的解为：\n");

X[n-1]=a[n-1][n]/a[n-1][n-1];

for(i=n-2;i>=0;i--)

{

s=a[i][n];

for(j=i+1;j<n;j++)

s=s-a[i][j]*X[j];

X[i]=s/a[i][i];

}

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

for(j=0;j<n+1;j++)

{

printf("%12f",a[u][j]);

c++;

if(c%(n+1)==0)

printf("\n");

}

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

{

printf("x(%d)=%.4f\n",i+1,X[i]);

if(i==n-1)

printf("\n");

}

system("pause");

return

0;

}

2、运行程序，方程的解为：x1x2x3x4x5x6x7x8x9-0.289230.34544-0.71281-0.22061-0.430400.15431-0.0578230.201050.290233、感想：首先，通过本上机题，提高了我运用循环语句的能力；其次，是我对各种排序算法有了更深的理解；最后，是我对二维数组有了更好的把握，并且认识到MATLAB在处理矩阵运算上的优势。在上机的过程中，加深了对Gauss消去法的认识。4、附图：图1输入增广矩阵习题437．〔上机题〕3次样条插值函数〔1〕编制求第一型3次样条插值函数的通用程序；〔2〕汽车门曲线型值点的数据如下：i012345678910Xi012345678910Yi2.513.304.044.705.225.545.785.405.575.705.80端点条件为y0'=0.8，y10'=0.2，用所编程序求车门的答：〔1〕通用程序：#

include<stdio.h>

#

include<math.h>

#

include<stdlib.h>

#

define

delta

1e-6

#

define

N

100

int

main()

{

int

i,j,t,r,n,u,c=0;

float

p,L,max,s,h,D1,D2;

float

x[N],M[N],y[N];

float

a[N][N+1]={0};

/*输入插值条件*/

printf("请输入插值节点的个数：\n");

scanf("%d",&n);

printf("请从左到右输入插值节点：\n");

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

{

scanf("%f",&x[i]);

}

printf("请输入插值节点对应的函数值：\n");

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

{

scanf("%f",&y[i]);

}

printf("请输入插值初始条件：\n");

printf("y'(0)=");

scanf("%f",&D1);

printf("y'(n)=");

scanf("%f",&D2);

/*产生弯矩方程组增广矩阵系数*/

a[0][1]=1;

a[n-1][n-2]=1;

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

a[i][i]=2;

for(i=1;i<n-1;i++)

{

h=(x[i]-x[i-1])/(x[i+1]-x[i-1]);

a[i][i-1]=h;

a[i][i+1]=1-h;

}

a[0][n]=6*((y[n-2]-y[n-1])/(x[n-2]-x[n-1])-D2)/(x[n-2]-x[n-1]);

a[n-1][n]=6*((y[1]-y[0])/(x[1]-x[0])-D1)/(x[1]-x[0]);

for(i=1;i<n-1;i++)

{

a[i][n]=6*((y[i+1]-y[i])/(x[i+1]-x[i])-(y[i]-y[i-1])/(x[i]-x[i-1]))/(x[i+1]-x[i-1]);

}

printf("请确认弯矩方程的增广矩阵为：\n");

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

{

for(j=0;j<n+1;j++)

printf("%3.2f\t",a[i][j]);

printf("\n");

}

for(j=0;j<n-1;j++)

{

{

max=fabs(a[j][j]);

r=j;

}

for(i=j+1;i<n;i++)

if(fabs(a[i][j])>max)

{

max=a[i][j];

r=i;

}

if(max<delta)

printf("矩阵奇异\n");

for(t=j;t<n+1;t++)

{

p=a[j][t];

a[j][t]=a[r][t];

a[r][t]=p;

}

for(i=j+1;i<N;i++)

{

L=a[i][j]/a[j][j];

for(t=j;t<n+1;t++)

a[i][t]=a[i][t]-L*a[j][t];

}

}

/*输出弯矩方程的计算结果*/

printf("整理后的增广矩阵为：\n");

M[n-1]=a[n-1][n]/a[n-1][n-1];

for(i=n-2;i>=0;i--)

{

s=a[i][n];

for(j=i+1;j<n;j++)

s=s-a[i][j]*M[j];

M[i]=s/a[i][i];

}

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

for(j=0;j<n+1;j++)

{

printf("%3.2f\t",a[u][j]);

c++;

if(c%(n+1)==0)

printf("\n");

}

printf("弯矩方程解向量为：\n");

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

{

printf("M(%d)=%f\n",i+1,M[i]);

if(i==n-1)

printf("\n");

}

/*计算插值函数*/

printf("样条函数如下：\n");

for(i=0;i<n-1;i++)

{

printf("S(x)=%f+%f(x-%.2f)+%f(x-%.2f)^2+%f(x-%.2f)^3\t%2.1f<x<%2.1f",y[i],((y[i+1]-y[i])/(x[i+1]-x[i])-(M[i]/3+M[i+1]/6)*(x[i+1]-x[i])),x[i],M[i]/2,x[i],(M[i+1]-M[i])/6/(x[i+1]-x[i]),x[i],x[i],x[i+1]);

printf("\n");

}

printf("S(i+0.5)的值如下：\n");

for(i=0;i<n-1;i++)

{

printf("S(%d+0.5)=%f",i,y[i]+((y[i+1]-y[i])/(x[i+1]-x[i])-(M[i]/3+M[i+1]/6)*(x[i+1]-x[i]))*(i+0.5-x[i])+M[i]/2*(i+0.5-x[i])*(i+0.5-x[i])+(M[i+1]-M[i])/6/(x[i+1]-x[i])*(i+0.5-x[i])*(i+0.5-x[i])*(i+0.5-x[i]));

printf("\n");

}

system("pause");

return

0;

}

〔2〕打印出3次样条插值函数S(x),并打印出S(i+0.5)，〔i=0，1，…，9〕1、样条函数如下：S2、S(i+0.5)的值如下：Xi0.51.52.53.54.55.56.57.58.59.5Yi2.89113.68314.38024.98855.38325.72385.59415.43115.65515.7497〔3〕附图：输入插值条件生成弯矩方程的增广矩阵输出插值函数及计算结果习题623．〔上机题〕常微分方程初值问题数值解常微分方程初值问题数值解〔1〕编制RK4方法的通用程序；〔2〕编制AB4方法的通用程序〔由RK4提供初值〕；〔3〕编制AB4-AM4预测校正方法通用程序〔由RK4提供初值〕；〔4〕编制带改良的AB4-AM4预测校正方法通用程序〔由RK4提供初值〕；〔5〕对于初值问题y取步长h=0.1，应用〔1〕~〔4〕中的四种方法进行计算，并将计算结果和精确解yx=3/(1+x〔6〕通过本上机题，你能得到哪些结论？答：〔1〕RK4通用程序：#include<stdio.h>

#include<stdlib.h>

#include<math.h>

#define

H

0.1

int

main()

{

float

x,y,r=3;

float

k1,k2,k3,k4;

float

fun1(float,float);

int

i,n;

printf("请输入x,y的初值：\n");

scanf("%f%f",&x,&y);

printf("请输入运算次数：\n");

scanf("%d",&n);

printf("运算结果为：\n");

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

{

printf("x%-2d=%.2f,\ty%-2d=%f,\ty(%-2d)=%f,\tr(%-2d)=%f",i,x,i,y,i,r,i,r-y);

printf("\n");

k1=fun1(x,y);

k2=fun1(x+H/2,y+H*k1/2);

k3=fun1(x+H/2,y+H*k2/2);

k4=fun1(x+H,y+H*k3);

x=x+H;

y=y+H*(k1+2*k2+2*k3+k4)/6;

r=3/(1+x*x*x);

}

system("pause");

return

0;

}

float

fun1(float

x,float

y)

{

float

f;

f=-x*x*y*y;

return

f;

}

〔2〕AB4通用程序：#include<stdio.h>

#include<stdlib.h>

#include<math.h>

#define

H

0.1

int

main()

{

float

x,y,r=3;

float

k1,k2,k3,k4,t[100]={0};

float

fun1(float,float);

int

i,j,n;

printf("请输入x,y的初值：\n");

scanf("%f%f",&x,&y);

printf("请输入运算次数：\n");

scanf("%d",&n);

printf("运算结果为：\n");

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

{

if(i<5)

{

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

{

printf("x%-2d=%.2f,\ty%-2d=%f,\ty(%-2d)=%f",j,x,j,y,j,r);

printf("\n");

k1=fun1(x,y);

k2=fun1(x+H/2,y+H*k1/2);

k3=fun1(x+H/2,y+H*k2/2);

k4=fun1(x+H,y+H*k3);

t[j]=y;

x=x+H;

y=y+H*(k1+2*k2+2*k3+k4)/6;

r=3/(1+x*x*x);

}

x=x-H;

}

t[i]=t[i-1]+H*(55*fun1(x,t[i-1])-59*fun1(x-H,t[i-2])+37*fun1(x-2*H,t[i-3])-9*fun1(x-3*H,t[i-4]))/24;

x=x+H;

r=3/(1+x*x*x);

printf("x%-2d=%.2f,\ty%-2d=%f,\ty(%-2d)=%f,\tr(%-2d)=%f",i,x,i,t[i],i,r,i,r-t[i]);

printf("\n");

}

system("pause");

return

0;

}

float

fun1(float

x,float

y)

{

float

f;

f=-x*x*y*y;

return

f;

}

〔3〕AB4-AM4通用程序：#include<stdio.h>

#include<stdlib.h>

#include<math.h>

#define

H

0.1

int

main()

{

float

x,y,r=3;

float

k1,k2,k3,k4,t[100]={0},Y[100]={0};

float

fun1(float,float);

int

i,j,n;

printf("请输入x,y的初值：\n");

scanf("%f%f",&x,&y);

printf("请输入运算次数：\n");

scanf("%d",&n);

printf("运算结果为：\n");

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

{

if(i<5)

{

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

{

printf("x%-2d=%.2f,\ty%-2d=%f,\ty(%-2d)=%f",j,x,j,y,j,r);

printf("\n");

k1=fun1(x,y);

k2=fun1(x+H/2,y+H*k1/2);

k3=fun1(x+H/2,y+H*k2/2);

k4=fun1(x+H,y+H*k3);

Y[j]=y;

x=x+H;

y=y+H*(k1+2*k2+2*k3+k4)/6;

r=3/(1+x*x*x);

}

x=x-H;

}

t[i]=Y[i-1]+H*(55*fun1(x,Y[i-1])-59*fun1(x-H,Y[i-2])+37*fun1(x-2*H,Y[i-3])-9*fun1(x-3*H,Y[i-4]))/24;

Y[i]=Y[i-1]+H*(9*fun1(x+H,t[i])+19*fun1(x,Y[i-1])-5*fun1(x-H,Y[i-2])+fun1(x-2*H,Y[i-3]))/24;

x=x+H;

r=3/(1+x*x*x);

printf("x%-2d=%.2f,\ty%-2d=%f,\ty(%-2d)=%f,\tr(%-2d)=%f",i,x,i,Y[i],i,r,i,r-Y[i]);

printf("\n");

}

system("pause");

return

0;

}

float

fun1(float

x,float

y)

{

float

f;

f=-x*x*y*y;

return

f;

}

〔4〕改良的AB4-AM4通用程序：#include<stdio.h>

#include<stdlib.h>

#include<math.h>

#define

H

0.1

int

main()

{

float

x,y,r=3;

float

k1,k2,k3,k4,t[100]={0},Y[100]={0},R[100]={0};

float

fun1(floa