東南大學(xué)數(shù)值分析上機(jī)作業(yè)

習(xí)題117．〔上機(jī)題〕舍入誤差與有效數(shù)設(shè)，其精確值為?！?〕編制按從大到小的順序，計(jì)算的通用程序?！?〕編制按從小到大的順序，計(jì)算的通用程序。〔3〕按兩種順序分別計(jì)算，，，并指出有效位數(shù)?！簿幹瞥绦驎r(shí)用單精度〕〔4〕通過本上機(jī)題，你明白了什么？答：1、從大到小順序：#include

<stdio.h>

int

main()

{

float

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

printf("請(qǐng)輸入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("請(qǐng)輸入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、計(jì)算結(jié)果：N101010從大到小S0.7400490.7498520.749852從小到大S0.7400500.7499000.749999精確值0.7400490.74990.749999有效位數(shù)從大到小S643從小到大S5464、感想：算法的不同導(dǎo)致計(jì)算結(jié)果的誤差不同，一定要好好設(shè)計(jì)算法，盡量使算法精確。通過上面的例子，可以看出按不同的計(jì)算順序所得到的結(jié)果是不同的，按從大到小的順序計(jì)算的值與精確值有較大的誤差，而按從小到大的順序計(jì)算的值與精確值吻合。計(jì)算機(jī)在進(jìn)行數(shù)值計(jì)算時(shí)會(huì)出現(xiàn)“大數(shù)吃小數(shù)〞的現(xiàn)象，導(dǎo)致計(jì)算結(jié)果的精度有所降低，使用計(jì)算機(jī)進(jìn)行同號(hào)數(shù)的加法時(shí)，采用絕對(duì)值較小者先加的算法，得到的結(jié)果的相對(duì)誤差較小。5、附圖：圖SEQ圖\*ARABIC1從大到小算法及結(jié)果圖SEQ圖\*ARABIC2從小到大算法及結(jié)果習(xí)題220．〔上機(jī)題〕Newton迭代法〔1〕給定初值及容許誤差，編制Newton法解方程根的通用程序?！?〕給定方程，易知其有三個(gè)根，，。1．由Newton方法的局部收斂性可知存在，當(dāng)時(shí)，Newton迭代序列收斂于根。試確定盡可能大的。2．試取假設(shè)干初始值，觀察當(dāng)，，，，時(shí)Newton序列是否收斂以及收斂于哪一個(gè)根?！?〕通過本上機(jī)題，你明白了什么？答：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("請(qǐng)輸入初值：");

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、確定盡可能大的：〔由于函數(shù)是奇函數(shù)，故僅需在正半軸求〕#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;

}

運(yùn)算結(jié)果為：0.7745963、測(cè)試是否在給定區(qū)間收斂：在五個(gè)區(qū)間內(nèi)分別取初值-100000，-0.8，0.5，0.9，110000測(cè)試是否收斂及收斂到哪個(gè)根，結(jié)果如下表：初值-100000-0.80.50.9110000是否收斂是是是是是收斂于-3-30334、感想：選取的迭代初值不同，則迭代序列將有可能收斂于不同的根，牛頓法對(duì)迭代初值的要求較高。有些區(qū)間收斂范圍相當(dāng)大，有些卻非常小。另外，我發(fā)現(xiàn)容許誤差的大小對(duì)迭代結(jié)果的精度有很大的影響，如果想得到比擬精確的結(jié)果，要將容許誤差設(shè)置小一點(diǎn)。5、附圖：圖1利用牛頓法求方程的近似解圖2確定收斂范圍習(xí)題339．〔上機(jī)題〕列主元Gauss消去法對(duì)于某電路的分析，歸結(jié)為求解線性方程組RI=V，其中R為9階矩陣，參見課本127頁(yè)VT=〔1〕編制解n階線性方程組Ax=b的列主元三角分解法的通用程序；〔2〕用所編制的程序解線性方程組RI=V，并打印出解向量，保存五位有效數(shù)；〔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("請(qǐng)輸入系數(shù)矩陣的階數(shù)：\n");

scanf("%d",&n);

printf("輸入的增廣矩陣系數(shù),中間用空格隔開：\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、運(yùn)行程序，方程的解為：x1x2x3x4x5x6x7x8x9-0.289230.34544-0.71281-0.22061-0.430400.15431-0.0578230.201050.290233、感想：首先，通過本上機(jī)題，提高了我運(yùn)用循環(huán)語(yǔ)句的能力；其次，是我對(duì)各種排序算法有了更深的理解；最后，是我對(duì)二維數(shù)組有了更好的把握，并且認(rèn)識(shí)到MATLAB在處理矩陣運(yùn)算上的優(yōu)勢(shì)。在上機(jī)的過程中，加深了對(duì)Gauss消去法的認(rèn)識(shí)。4、附圖：圖1輸入增廣矩陣習(xí)題437．〔上機(jī)題〕3次樣條插值函數(shù)〔1〕編制求第一型3次樣條插值函數(shù)的通用程序；〔2〕汽車門曲線型值點(diǎn)的數(shù)據(jù)如下：i012345678910Xi012345678910Yi2.513.304.044.705.225.545.785.405.575.705.80端點(diǎn)條件為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("請(qǐng)輸入插值節(jié)點(diǎn)的個(gè)數(shù)：\n");

scanf("%d",&n);

printf("請(qǐng)從左到右輸入插值節(jié)點(diǎn)：\n");

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

{

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

}

printf("請(qǐng)輸入插值節(jié)點(diǎn)對(duì)應(yīng)的函數(shù)值：\n");

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

{

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

}

printf("請(qǐng)輸入插值初始條件：\n");

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

scanf("%f",&D1);

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

scanf("%f",&D2);

/*產(chǎn)生彎矩方程組增廣矩陣系數(shù)*/

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("請(qǐng)確認(rèn)彎矩方程的增廣矩陣為：\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];

}

}

/*輸出彎矩方程的計(jì)算結(jié)果*/

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

}

/*計(jì)算插值函數(shù)*/

printf("樣條函數(shù)如下：\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次樣條插值函數(shù)S(x),并打印出S(i+0.5)，〔i=0，1，…，9〕1、樣條函數(shù)如下：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〕附圖：輸入插值條件生成彎矩方程的增廣矩陣輸出插值函數(shù)及計(jì)算結(jié)果習(xí)題623．〔上機(jī)題〕常微分方程初值問題數(shù)值解常微分方程初值問題數(shù)值解〔1〕編制RK4方法的通用程序；〔2〕編制AB4方法的通用程序〔由RK4提供初值〕；〔3〕編制AB4-AM4預(yù)測(cè)校正方法通用程序〔由RK4提供初值〕；〔4〕編制帶改良的AB4-AM4預(yù)測(cè)校正方法通用程序〔由RK4提供初值〕；〔5〕對(duì)于初值問題y取步長(zhǎng)h=0.1，應(yīng)用〔1〕~〔4〕中的四種方法進(jìn)行計(jì)算，并將計(jì)算結(jié)果和精確解yx=3/(1+x〔6〕通過本上機(jī)題，你能得到哪些結(jié)論？答：〔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("請(qǐng)輸入x,y的初值：\n");

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

printf("請(qǐng)輸入運(yùn)算次數(shù)：\n");

scanf("%d",&n);

printf("運(yùn)算結(jié)果為：\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("請(qǐng)輸入x,y的初值：\n");

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

printf("請(qǐng)輸入運(yùn)算次數(shù)：\n");

scanf("%d",&n);

printf("運(yùn)算結(jié)果為：\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("請(qǐng)輸入x,y的初值：\n");

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

printf("請(qǐng)輸入運(yùn)算次數(shù)：\n");

scanf("%d",&n);

printf("運(yùn)算結(jié)果為：\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