浙江大学数值分析C语言编程习题

C语言编程习题

第二章

习题2・2

5.用二分法编程求6x4-40x2+9=0的所有实根。

include<stdio.h>

#_nclude<math.h>

defineN10000

doubleA,B,C;

doub1ef(doublex)

|

return(A*x*x*x*x+B*x*x+C);

)

voidBM(doublea,doubleb,doubleepsl,doubleeps2)

(

inik;

doublex,xe;

doublevaluea=f(a);

doublevalueb=f(b);

if(valuea>0&&valueb>0||valuea<0&&valueb<0)return;

prinlf("Findingrootintherange:[%.31f,%.31f]\n”,a,b);

for(k=l;k<=N;k++){

x=(a+b)/2;

xe=(b-a)/2;

if(fabs(xe)<eps2)fabs(f(x)Xepsl){

printf(/zThexvalueis:%g\n',x);

printf("f(x)=%g\n\n”,f(x));

return;

)

if(f(a)*f(x)<0)b=x;

elsea=x;

)

printf(?，Noconvergence!\n*);

)

intmain()

(

doublea,b,epsl,eps2,step,start;

printf("PleaseinputA,B,C:\n*);

scanfC%lf%lf&B,&C);

printf("Pleaseinputa,b,step,epsl,eps2:\n");

scant-C%lf%lf%lf%lf%ir,&a,&b,&step,&epsl,&eps2);

for(slart=a;(start+step)<=b;start+=step){

doubleleft=start;

dcMihlpright=start+atop：

BM(left,right,epsl,eps2);

return0;

)

运行：

PleaseinputA,B,C:

6-409

Pleaseinputa,b,step,epsl,eps2:

-10101le-5le-5

Findingrootintherange:[-3.000,-2.000]

Thexvalueis:-2.53643

f(x)=-0.00124902

Findingrootintherange:[-1.000,0.000]

Thexvalueis:-0.482857

f(x)=0.00012967

Findingrootintherange:[0.000,1.000]

Thexvalueis:0.482857

f(x)=0.00012967

Findingrootintherange:[2.000,3.000]

Thexvalueis:2.53643

f(x)=-0.00124902

有时假设把判别语句

if(fabs(xe)<eps2Ifabs(f(x))<epsl)

改为

if(fabs(xe)<cps2&&fabs(f(x)Xcpsl)

会提高精度，对同一题运行结果：

Findingrootintherange:[-3.000,-2.000]

Thexvalueis:-2.53644

f(x)=-4.26496e-007

Findingrootintherange:[-1.000,0.000]

Thexvalueis:-0.482861

f(x)=-7.3797e-006

Findingrootintherange:[0.000,1.000]

Thexvalueis:0.482861

f(x)=-7.3797e-006

Findingrootintherange:[2.000,3.000]

Thexvalueis:2.53644

f(x)=-4.26496e-007

习题2・3

5.请用埃特金方法编程求出x=tgx在4.5(弧度)附近的根。

#include<stdio.h>

ttmclude<math.h>

^defineN100

ticlpfinpPT3.1415926

voidSM(doublexO,doubleeps)

(

intk;

doublex;

doublexl,x2;

for(k=l;k<=N;k++){

xl=sin(xO)/cos(xO);

x2=sin(xl)/cos(xl);

x=xO-(xl-xO)*(x1-xO)/(x2-2*x1+xO);

if(fabs(x-xO)<eps){

printf(''Thexvalueis:%g\n',x);

return;

)

xO=x;

)

printf(//Noconvegence!\n*);

1

intmainO

(

doubleops,xO;

printf("Pleaseinputeps,xO:\n");

scant'Cz%lf&eps,&x0);

SM(xO,eps);

return0;

)

运行：

Pleaseinputeps,xO:

le-54.5

Thexvalueis:4.49341

习题2.4

11.请编出用牛顿法求复根的程序，并求出

P(z)=z4-3z3+20z2+44z+54=0

接近于zo=2.5+4.5i的零点。

#include<cstdio>

#inchide<cmath>

#defineMAX_TIMES1(X)()

typedefstruct{

doublereal,image;

}COMPLEX;

COMPLEXAa[5]={{54,0},{44,0),{20,0},{-3,0},{1,0}};

returnsqrt(a.real*a.real+a.image*a.image);

)

COMPLEX((COMPLEXx)

{

in(i;

COMPLEXresult=zero;

for(i=0;i<5;i++){

result=add(result,multi(Aa[i],times(x,i)));

)

returnresult;

)

COMPLEXff(COMPLEXx){

inti;

COMPLEXresuh=zero;

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

rcsult=add(rcsult,niulti(BbliJ,timcs(x,i)));

)

returnresult;

)

intmain(){

COMPLEXzO,zl,result;

doublex,y;

intk;

printf("pleaseinputx.y\n");

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

zO.real=x;zO.image=y;

for(k=0;k<MAX_TIMES;k++){

zl=sub(ract(zO,Div(f(zO),ff(zO)));

result=f(zl);

if(distance(z(),zl)<epsl||complex_abs(result)<eps2){

printf("Therootis:z=(%.3f+%.3fi),f(z)=(%e+%ei)\n",zl.real,zl.image,result.real,

result,image);

return0;

}

zO=zI;

I

printf("Noconvergence!\n");

return0;

)

运行：

pleaseinputx,y

2.54.5

Therootis:

z=(2.471+4.64li),f(z)=(-1.122705e-007+-1.910245e-007i)

习题2・6

2.请编程用劈因子法求高次方程d+/+5/+4;v+4=0的所有复根。

#include<cstdio>

#inchide<cmath>

intmain()

{

floatb[20J,c[20];

floatdelta;

floatcps;

intprint=O;

floatfruO,frul,sO,s1,frvO,frv1;

floatdeltau,deltav;

floatu,v;

floata[20J;

intn,j;

floatrO,rl;

floatr;

inti,k;

printf("PleaseinputmaxexpC'nentAn");

scanf("%d".&n);

printf("Pleaseinputcpslion:\n");

scanf("%f',&eps);

printf("Pleaseinputthecoefficient:\n");

for(i=0;i<=n;i++)scanf("%f,",&a[i]);

for(j=0;j<=n;j++)printf(1,+%.3fXA%d",aLi],n-j);

printf("\n");

for(k=l;k<=2;k++){

u=a[n-lJ/aln-2];

v=a[n]/a[n-2];

if(k==2){

u=0;v=4;

I

while(l){

b[0]=a[0];

b[l]=all]-u*b[0];

forG=2;j<=n;j++){

bU]=aU]-u*bU-l]-v*b[j-2];

)

rO=b[n-l];

rl=b[n]+u*b[n-l];

c[0]=b[0];

c[lJ=bll]-u*b[OJ;

forG=2;j<=n-2;j++)

c[j]=b|j]-u*c[j-1]-v*c[j-2];

s0=c[n-3];

sl=c[n-2]+u*c[n-3];

fruO=u*sO-sl;

frul=v*sO;

frvO=-sO;

frv1=-s1;

deltau=-(frvI*r()-frvO*rl)/(fruO*frv1-fru1*frvO);

del(av=-(frul^rO-fruO*riyCfrvO*frul-frv1*fruO);

u=u+dcitau;

v+=deltav;

delta=u*u-4*v;

if((fabs(dcltau)<cp;s)&&(fabs(dcUav)<cps))

{

print=l;

printf("foundroots:");

r=-u/2;

i=(int)sqrl(fabs(deha))/2;

if(dclta<0){

printf("%,3f+%.3fi\t",r,fabs((double)i));

printf("%,3f-%.3fi\n",r,fabs((double)i));

)

else(

printf("%,3At\tu,r+i)；

printf("%.3An",r-i);

I

if(n==l)printf("Foundrooi:%.3f\n';a[l]/a[0]);

if(k==2)return0;

)

if(k==l&&print)break;

}/*endwhile*/

}/*endfor*/

return0;

)

运行：

Pleaseinputmaxexponent:

4

Pleaseinputepslion:

le-5

Pleaseinputthecoefficient:

11544

+L000X人4+L000X人3+5.000XA2+4.000XA1+4.000XA0

习题34

5.请编程用列全主元高斯一约当消去法求矩阵A,B的逆矩阵。

#include<cstdio>

#include<cmath>

#defineMAX10

intROW;

staticfloatA[MAX+1][2*MAX+1];

voidputout()

inti,j;

printf("\nThereversematrixis:\n");

for(i=l;i<=ROW;i++)|

fbr(j=ROW+l;jv=2*ROW;j++){

printf(H%.3fH,A[i][j]);

}

printf("\n");

I

)

intget()

{

floatmax[MAX+l];

intcount_r[MAX+1];

floattmp[2*MAX+l];

floatpassI,pass2;

inti,j,k;

for(i=l;i<=ROW;i++){

max[i]=A[i][i];

count_rlij=i;

for(k=i;k<=ROW;k++){

if(max[i]<A[k][i]){

max[i]=A[k][iJ;

count_r[i]=k;

)

}/***********findthemaxone**********/

if(max[i]==0)return-1;

if(count_r[i]!=i){

fdr(k=1；k<=2*ROW;k++)

tmp[k]=A[count_i[i]][k];

A[countjr[i]][k]=A[i][k];

A[i][k]=tmp[k];

I

}/*********chantherows*************/

passl=A[i][ij;

fbr(k=l;k<=2*ROW;k++){

A[i][k]/=passl;

for(j=l;j<=ROW;j++){

if(j==i)continue;

else{

pass2=A[j][i]；

for(k=l;k<=2:}ROW：k++){

AU]M=AU][k]-pass2*A[i][k];

।/************getsimple**************/

putoutO;

return0;

intmain()

{

inii,j;

floattmp=-1;

floatp;

printf("pleaseinputthenumberofROW:");

scanf("%d",&ROW);

printf("\n");

if(ROW>=MAX){

primf("ihenumberofrowistoobig!\n");

return0;

)

printf("PleaseinputthematrixA:\n");

for(i=l;i<=ROW;i++){

for(j=l；j<=ROW;j++){

scanf("%f；\&p);

A[i]|j]=p；

if(tmp<fabs(A[i][j]))tmp=fabs(A[i]|j]);

)

}

if(trnp==O){

printf("Noinverse!");

return0;

I

for(i=l;i<=ROW;i++){

AliJ[i+ROW]=l;

}

returnget();

)

运行：

pleaseinputthenumberofROW:

3

PleaseinputthematrixA:

-385

2-74

19-6

Thereversematrixis:

0.0260.3960.285

0.0680.0550.094

0.1060.1490.021

pleaseinputthenumberofROW:

4

PleaseinputthematrixA:

21-3-1

3107

-124-2

10-15

Thereversematrixis:

-0.0470.588-0.271-0.941

0.388-0.3530.4820.765

-0.2240.294-0.035-0.471

-0.035-0.0590.0470.294

习题3.2

4.编程用追赶法解以下三对角线方程组。

#include<stdio.h>

#includc<math.h>

#defineROW5

staticfloatA|ROW+1]={0,0,0,0,0.0};

sialicGoatB[ROW+1]={0,0,0,0,0.0};

staticfloatC[ROW+l]={0A0A0.0);

staticfloatX|ROW+IJ={0,0,0,0,0,0);

intmain()

(

inti;

printf("PlcascinputtheA:\n");

for(i=2;i<=ROW;i++)scanfC'%f,u,&A[i]);

printf("PleaseinputtheB:\n");

for(i=l;i<=ROW;i++)scanf("%f,",&B|iJ);

printf("PleaseinputtheC:\n");

for(i=l;i<=R0W-l;i++)scanf(”％f,”,&C[i])；

printf("PleaseinputtheF:\n");

for(i=l;i<=R0W;i++)scanf(”％f,”,&X[i]);

C[I]=C[I]/B[1];

for(i=2;i<=ROW-l;i++)

C[iJ=C[i]/(Blij-A[i]*C[i-lJ);

X[1]=X[1]/B[1];

for(i=2;i<=ROW;i++)

X[i]=(X[i]-A[i]*X[i-l])/(B[i]-A[i]*C[i-l]);

for(i=ROW-l;i>=l;i-)

X[i]=X[i]-C[i]*X[i+l];

printf("Thevalueis:\n");

for(i=l;i<=ROW;i++)printf(nx%d=%.3f^',,i,X[i]);

return0;

)

运行：

PleaseinputtheA:

-1-1-1-1

PleaseinputtheB:

44444

PleaseinputtheC:

-1-1-1-1

PleaseinputtheF:

100200200200100

Thevalueis:

xl=46.154

x2=84.615

x3=92.308

x4=84.615

x5=46.154

习题4・3

1.用高斯一赛德尔方法编程解以下线性方程组，要求当

口供+1)伏)时迭代终止。

#include<cstdio>

#include<cmath>

floatx[6];

floaty[6];

floateps;

floata[6][6]={

{4,-1,0,-1,0,0},

{-l,4,-kO,-l,O),

[0,-1,4,0,0,-1),

{-1,0,0,4,-1,0),

{0,-1,0,-1,4,-1),

{0.0,-1,0,-1,4}

}；

floatb[6]={0,5,0,6.-2,6);

intgs()

inti,k,j;

floats;

for(i=0;i<6;i++)y[il=x[i]=0;

for(k=0;k<20;k++){

for(i=0;i<6;i++){

s=0;

forG=0;j<6：j++){

if(j!=i)s=s+a[i]|j]*y|j];

)

y[i]=(b[i]-s)/a[i][i];

I

for(i=0;(i<6)&&(abs(y[i]-x[i])<eps);i++);

if(i>=5)break;

for(j=0;j<6;j++)x[j]=yfj];

if(k>20){

printf("Noconvcrgcncc./n");

return-1;

I

}

returnI;

)

intmain()

{

inttag,i,m;

prinif("\mhemaxcircles：");

scanf("%d",&m);

printf("eps=");

scanf("%lf',&eps);

tag=gs();

if(tag>0){

for(i=0;i<=5;i++){

x[i]=y[i]；

printf("x(%d)=%.3e\n,,,i+1,x[i]);

)

}

return0;

)

运行：

nthemaxcircles=

10

eps=

le-5

x(l)=1.000e+000

x(2)=2.000e+000

x(3)=1.000c+000

x(4)=2.000e+000

x⑸=L000e+000

x(6)=2.000e+000

习题4-4

#include<cmath>

#includc<cstdio>

intagsdl(doublea[][101,doubleb[1.intn,

doublex(],doublccps.doublc\v,intm)

(

inti,j,pcount;

doublep,t,s,q;

for(i=0;i<10;i++)x[i]=0;

p=eps+1.0;

for(pcount=0;p>=eps&&pcount<m;pcount++){

p=0.0;

for(i=0;i<n;i«»){

t=x[i];s=0.0;

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

if(j!=i)s=s+a[i]|j]*x[j];

x[i]=(b[i]-s)/a[i][i];

q=w*fabs(xfi]-t);

if(q>p)p=q;

1

)

returnpcouni;

}

intniain()

{

inti,j,n,m,count=0;

doubleeps,w[3];

staticdoublea[10][10];

staticdoublex[10],b[10];

piintf("\nthemaxcircles=");

scanf("%d",&m);

printf("n=");

scanf("%d",&n);

printf("\nmatrixA:\n");

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

for(j=0;j<n;j++)scanf("%lf\&a[i][j]);

printf("\nvectorb:\n");

for(i=0;i<n;i++)scanf("%lf',&b[ij);

printf("eps=");

scanf("%lf',&eps);

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

printf("w(%d)=H,i);

scanf("%If\&w[i]);

I

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

printf("w(%d)=%.3f\n',,j,w[j]);

if(agsdl(a,b,n,x,eps,wLj],in)<in){

for(i=0;i<n;i++)printf("x(%d)=%.3e”,i,x[i]);

printfC'Xn");

)

else

{

printfCfaibn1');

return0;

)

运行：

themaxcircles=

10

n=

3

matrixA:

4-10

-14-1

0-14

vectorb:

14-3

eps=

le-5

w(0)=

1.03

w(l)=

1

w(2)=

1.1

vv(0)=1.030

x(0)=5.000e-001x(l)=1.000e+000x(2)=-5.000e-001

w⑴=1.000

x(0)=5.000e-001x(l)=1.000e+000x(2)=-5.000e-001

w(2)=1.100

x(0)=5.000e-001x(l)=1.000e+000x(2)=-5.000e-001

第5章

习题5・1

2.编程求以下矛盾方程组的解:

#include<fstream.h>

#include<stdio.h>

#include<stdlib.h>

doubleeps=1.0e-5;//errorrange

double**a;//savematrixA

double**ab;//saveaugmentedmatrixAB

double**ta;//savematrixATA

double*b：//savevectorsB

double*tb;//saveATB

double**ca;//savematrixA

double*cb;//savevectorsB

double*x;//saveresultvectorsX

int*fr;//rowflag

int*fc^/columnflag

intm,cm;//rownumofmatrixA

intn;//columnnumofmatrixA

intr;

voidinit();

iniG_J();

intd_O(int);

voidchange();

voidg_a_m();

voidshow();

voidshow1();

voidmain()

{iniI；

init();

t=l；

while(l)

{r=GJ();

if(m==r)break;

if(d_O(r)==O)break;

changeO;

t=0;

)

if(n==r)

{if(t)

{couivv”该线性方程组是非矛盾线性方程组，有唯一解："vvendl;

show。；

1

else

{cout<<”该线性方程组是矛盾线性方程组，有唯一最小二乘解："<<endl;

show。；

else

{if(t)

{coutvv”该线性方程组是非矛盾线性方程组，有普遍解：“vvendl;

show1();

)

else

{COUIV<”该线性方程组是矛盾线性方程组，有普遍最小二乘解："<<endl;

show1();

}

)

)

voidinit()

{//open(hedatafile-ab.txt

ifstreanifin("ab.txt");

if(!fin)

{cout«"Can'topenfileab.txttoread!"«endl;

exit(-l);

I

//readtherowandcolumninformationofmatrixA

fin»m»n;

cm=m;

//requireroomofmatrixA

a=(double**)malloc((m4-l)*sizeof(double*));

inti,j;

for(i=l;i<=m;i++)

a[i]=(doublc*)malloc((r+1)*sizcof(doublc));

//readmatrixA'sdatafromab.ixi

for(i=l;i<=m;i++)

for(j=l;j<=n;j++)

fin»a[i][j];

//initmatrixB

b=(double*)malloc((m+1)*sizeof(double));

for(i=l;i<=m;i++)

fin»b[i];

//closethedatafile

fin.close();

//backupmatrixAandB

ca=a;cb=b;

return;

}

//generateaugmentedmatrix

voidg_a_m()

{inti,j；

//requireroom

ab=(double**)malloc((m+1)*sizeof(double*));

for(i=l;i<=m;i++)

ab[i]=(double*)malloc((n+2)*sizeof(double));

//loaddata

for(i=l;i<=m;i++)

{for(j=l;j<=n;j++)

ab[i][j]=a[i][j];

ab[i][n+lj=b[i];

//inittheflagarrayfr(flagrow)andfc(flagcolumn)

fr=(int*)nialloc((m+l)*sizeof(int));

fc=(int*)malloc((n+l)*sizeof(int));

for(i=l;i<=m;i++)fr[i]=i;

for(j=l;j<=n;j++)fc|j]=j;

return;

}

//Gauss-Jordoneliminationmethod

intG_J()

{intk,found,min;

inti,j,pi,pj,ft；

doubletl,t2;

//generateaugmentedmatrix

g_a_m();

k=l;min=(m<n)?m:n;

//lookformainelement

whilc(k<=min)

{(i=found=pi=pj=0;

fbr(i=k;i<=m;i++)

for(j=k；j<=n;j++)

{//t2=abs(ab[fr[i]][fc[j]]);

t2=ab[fr[i]][fc|j]];

t2=(t2>O)?t2:(-t2);

if((t2>tl)&&(t2>eps))

{found=l;tl=t2;

Pi=i；Pj=j；

)

}//endofj

//endofi

if(!found)

returnk-l;

//foundamainelement

//fixtheflagarray-exchangetherowandcolumn

ft=fr[pi];fr[pi]=fr[k];fr[k]=ft;

ft=fc[pj];fc[pj]=fc[k];fc[k]=ft;

//rowstandardize

for(j=l;j<=n+l;j++)

ab[fr[kJJU]/=tl；

//elimination

for(i=l;i<=m;i++)

{if(i==k)continue;

t2=ab[fr[i]][fc[k]];

for(j=I;j<=n+l;j++)

ab[fr[i]]rj]-=ab[fr[k]][j]*t2;

)

//lookfornextmainelement

k++;

}//endofwhile

returnk-1;

)

//a=aTa,b=aTb

voidchange()

(

//requireroomforaTa

ta=(double**)malloc((n+1)*sizeof(double*));

inti;

for(i=l;i<=n;i।i)

tali]=(double*)nialloc((n+l)*sizeof(double));

//computingaTa

intj,k;

doublesum;

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

for(j=l;j<=n;j++)

{sum=0;

for(k=l;k<=m;k++i

sum+=a[k][i]*a[k][j];

ta[i][j]=sum;

)

//requirerooforaTb

tb=(double*)malloc((n+1)*sizeof(double));

//computingaTb

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

{sum=0;

for(k=l;k<=m;k++)

sum+=a[k][i]*b[k];

tb[i]=sum;

)

//freea,b

//a=aTa,b=aTb

a=ta;b=tb;

m二n;

return;

}

//iffromd(r+1)tod(m)areallzero,return0

//otherwise,return1

intd_O(intr)

{

intj;

inttO;

floattl;

t0=0;

for(j=r+1;j<=m;j++)

{tl=ab[frU]][n+l];

if(t1>eps||t1<-eps)

tO=l;

)

returntO;

}

//showtheaugmentedmatrixAB

voidshow()

{inti,j;

//outputandstoretheresultX[0..n-l]

x=(doublc*)malloc(n*sizcof(doublc));

for(j=l;j<=n;j++)

fbr(i=l;i<=m;i++)

if(ab[i][j]==l)

{xU-l]=ab[i][n-l];

cout«"X["«j-l«,,]="«xU-l]«endl;

}

//outputtheresidualvectors

cout«"该方程组的剩余向量是：”<«ndl;

coui〈〈"r=(";

doublesum;

for(i=l;i<=cm;i++)

{sum=0;

for(j=l;j<=n;j++)

sum+=ca[i][j]*x[j-l];

cout«cb[i]-sum;

if(i!=cm)cout«",";

)

cout«")T"«endl;

return;

I

voidshow1()

{inti,j,k.fbund;

for(i=l;i<=r;i++)

{cout«',X[,,«fc[i]«H]="«ab[fr[i]][n+l];

for(j=l;j<=n-r;j++)

cout«,,-X["«j«,']*(，，«ab[fr[i]][fc[r+j]]«")"；

cout«endl;

)

for(i=r+l:i<=n;i++)

cout«"X["«fc[i]«"J=X["«fc[i]«,']"«endl;

coutvv”是否需要特解?(y/n),;

charch;

cin»ch;

if((ch=-n')||(ch==,N'))return;

couivv"请输入参数："v<endl;

//inputvariables

x=(double:f:)malloc((n+1)*sizeof(clouble));

for(i=r+l;i<=n;i++)

{cout«"X[M«fc[i-r]«"l=";

cin»x(i];

)

//computingspecialresult

for(i=l;i<=r;i++)

{x[fc[i]]=ab[fr[i]][n+l];

fbr(j=l;j<=n-r;j++)

x[fc[i]l-=x[rrj]*ab[fr[i]][fc[rij]];

}

//outputspecialresult

coutvv”该方程组的特解是："«endl;

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

cout«"X[,,«i«,,]="«x[i]«endl;

//outputtheresidualvectors

couivv”该方程组的剩余向量是："vvendl;

cout«"r=(";

doublesum;

for(i=l;i<=cm;i++)

{sum=0;

for(j=l;j<=n;j++)

surn+=ca[ij[j]*x[j];

cout«cb[il-sum;

if(i!=cm)cout«'\";

}

cout«")T"«endl;

return;

}

运行：

ab.txt：

43

111

421

931

1641

4101826

该线性方程组是矛盾线性方程组，有唯一最小二乘解:

X[0]=0.5

X[l]=4.9

X[2]=-1.5

该方程组的剩余向量是：

r=(0.l,・0.3,0.3,-0・l)T

第六章

习题6・1

9.正弦函数表

Xk0.50.7().91.11.31.51.71.9

Sin冰0.47940.64420.78330.89120.96360.99750.99170.9463

编程用Lagrange插值多项式计算x()=0.6,0.8,1.0处的函数值sin(0.6),sin(0.8),

sin(l.O)的近似值"0.6)J(0.8)/(L0)。

#include<cstdio>

#include<cmath>

#include<cstdlib>

#defineN8

floatx[N]={0.5f,0.7f,0.9f,l.lf,1.3f,1.5f,1.7f,1.9f};

floatfTN]=(0.4794f,0.6442f,07833f,0.8912f,0.9636f,0.9975f,0.9917f,0.9463f);

floatlagrange(floatxO)

(

floats;

floatresult=0.0;

inlj,k;

for(k=0;k<N;k++){

s=1.0;

for(j=0;j<N;j++){

if(j!=k)s*=(xO-x[j])/(x[k]-x[j]);

1

result+=f[k]*s;

)

returnresult;

)

intmain(void)

(

floatxO,fO;

printf("InputthevalueofxO:");

scanf(H%r,&xO);

floatresult=lagrange(xO);

printf("xO=%.3ff(x0)=%.3fsin(x0)=%.3f\n",xO,result,sin(xO));

return0;

I

运行：

InputthevalueofxO:

x0=0.600f(x0)=0.565sin(x0)=0.565

InputthevalueofxO:

x0=0.800f(x0)=0.717sin(x0)=0.717

InputthevalueofxO:

x0=1.000f(x0)=0.841sin(x0)=0.841

习题6・5

7.某沿海货轮475水线的型值表,n=15

iXiV

11162.599.5

22950220.2355

35900589.83954

47675.8515950

588501260.2342

610881.2551900

7118802207.0194

813802.8992850

9147503124.1306

1017467.4813800

11177003853.2175

12206504390.3286

边界条件为

请编程序上机用三弯矩方程求三次样条函数第一类定解问题的解，

求xo=-9OO开始，每间隔AxulOO。,一共no=31个点上的函数值。单位毫米，当

插值点x落在[1162.5,29500]外面时,规定节点处的y值取大数10%

#include<cstdio>

#include<cmath>

#defineMAX10000000

#defineN15

doublex[]={ll62.5,2950,5900,7675.8515.8850,10881.255,11800.

13802.899,14750,17467.481,17700,20650,23600,26550,29500};

doubley[]={99.5,220.2355,589.83954,950,1260.2342,19(X),

2207.0194,2850,3124.1306,3800,3853.2175,4390.3286,4756.1508,

4976.2542,5028.7000};

intWINTH;

doubledy0=0.04770,dyn=0.(M)8810;

doubleM[N1;

doubleg[N],u[N],v[N],t[N];

doubleh[N];

voidcal(void)

inti;

for(i=0;i<N:i++)

v[O]/=t[O];

for(i=l;i<N-l;i++){

t[i]-=u[i-l]*v[i-l];

v[i]/=t[i]；

>

tlN-lJ-=ulN-2J*v[N-2J;

Mro]/=t[o]；

for(i=1;i<N;i++)M[i]=(M[i]-u[i-l]*M[i-l

for(i=N-2;i>=0;i-)M[i]-=v[i|*M[i+1];

)

intget_n(doubletemp)

{

inti;

if(temp<x[O]){

return0;

}elseif(temp>x[N-l]){

returnN;

}else{

for(i=l;i<N;i++)

if(temp>=x[i-1]&&temp<=x[i])

returni;

I

return0;

)

doublespline(doublek)

{

inti;

i=get_n(k);

if(i<I||i>=N)returnMAX;

else{

k=M[i-l]*(x[i]-k)*(x[i]-k)*(x[i]-k)/(6.0*h[i-l])

+M[i]*(k-x[i-1])*(k-x[i-l])*(k-x[i-1])/(6.0*h[i-1])

+(y[i-l]-M[i-l]*h[i-I]*h[i-1]/6.0)*(x[i]-k)/h[i-l]

+(y[i]-MLi]*h[i-1J*h[i-1J/6.0)*(k-x[i-1J)/h[i-1];

returnk;

)

)

voidoutput(void)

doublexO;

for(x0=-900;x()<=29500;x0+=WINTH){

printf("xO=%.21f,spline(x0)=%-16.2e\n",x0,spline(x0))；

)

)

voidinain()

inti;

prinlfC'pleaseinputdeltaX:\n");

scanf("%d'\&WINTH);

for(i=0;i<N;i++)t[i]=2;

for(i=0;i<N-1;i++)h[i]=x[i+1]-x[i];

for(i=0;i<N-2;i++)u[i]=h[i]/(h[i]+h[i+1]);

u[N-2]=1.0;

for(i=l;i<N-l;i++)v[i]=l-u[i-l];

v[O]=I.O;

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

g[i]=6/(h[i-I]+h[i])*((y[i-H]-y[i])/h[i]

g[0]=6.0*((y[1]-y[O])/h[O]-dyO)/h[O];

g[N-1]=6.0*(dyn-(y[N-1]-y[N-2])/h[N-2])/h[N-2];

cal();

output();

运行：

pleaseinputdeltaX:

1000

x0=-900.00,spline(x0)=1.00e+007

x0=1000.00,spline(x0)=1.00e+007

x0=1100.00,spline(x0)=1.00e+007

x0=2100.00,spline(x0)=1.54e+002

x0=3100.00,spline(x0)=2.34e+002

x0=4100.00,spline(x0)=3.36e+002

x0=5100.00,spline(x0)=4.65e+002

x0=6100.00,spline(x0)=6.24e+002

x0=7100.00,spline(x0)=8.19e+002

x0=8100.00,spline(x0)=1.06e+003

x0=9100.00,spline(x0)=1.33e+003

x0=10100.00,spline(x0)=1.64e+003

x0=11100.00,spline(x0)=1.97e+003

x0=12100.00,spline(x0)=2.31e+003

x0=13100.00,spline(x0)=2.63e+003

x0=14100.00,spline(x0)=2.94e+003

x0=15100.00,spline(x0)=3.22e+003

x0=16100.()0,spline(x0)=3.48e+003

x0=17100.00,spline(x0)=3.71e+003

x0=18100.00,spline(x())=3.94e+003

xO=19IOO.OO,spline(x0)=4.14e+003

x0=20100.00,spline(x0)=4.31e+003

x0=21100.00,spline(x0)=4.46e+003

x0=22100.00,spline(x0)=4.59e+003

x0=23100.00,spline(x0)=4.70e+003

x0=24100.00,spline(x0)=4.81e+003

x0=25100.00,spline(x0)=4.89e+003

x0=26100.00,spline(x0)=4.95e+003

x0=27100.00,spline(x0)=5.00e+003

x0=28100.00,spline(x0)=5.02e+003

x0=29100.00,spline(x0)=5.03e+003

第七章

习题7・2

2.对y=sinx在[0,90],取1度为间隔，分别用库函数sin(x)和用正交多项式对

sinx做四次曲线拟合，比拟拟合的结果与库函数的计算结果，并给出拟合曲线

的系数.

#include<cstdio>

#include<cstdlib>

#include<cmath>

#defineN91

#defmeNUM1()

#define03.1415926/180

#deHneK4

doubledelta(doubles[],intns,doublex[],doubley[Lintm);

doublePoly(doublep[],intn,doublexO)

(

registerintk;

doublef;

f=p(nl；

for(k=n-l;k>=0;k-){

f=f*xO+p[k];

}

returnf;

)

voidmulpoly(doublepl"intnl,doublep2[],intn2,doublep[],int*n)

(

intij；

*n=nl+n2;

for(i=0;i<=nl+n2;i++)p[i]=();

for(i=0;i<=nl;i++){

for(j=0y<=n2y++){

p[i+j]+=pl[i]*p2UJ；

)

}

)

voidaddpoly(cloublepl[],intnl,doublecl,doublep2[JJntn2,doublec2,doublep[],int*n)

inti;

intmin;

*n=(nl>=n2)?nl:n2;

min=(nl<=n2)?nl:n2;

for(i=0;i<=min;i++){

p[i]=pl[i]*cl+p2[i]*c2;

}

for(i=min+l;i<=nl;i++){

}

for(i=min+l;i<=n2;i++){

plij=p2[i]*c2;

}

)

voidzj(doublex[],doubley[],intm)

(

doublesum^unil,resu!t[N];

doublep[NUM][NUM],temppoly[NUM],a[NUM],temppolyl[2];

intn[NUMl,ntemp,nr;

inti,k;

doublealf[NUMl,beta[NUMl,temp4empl,eps=0.00001;

n[0]=0;

p[0][0]=l;

for(sum=(),i=0;i<in;i++){

sum+=y[i];

)

a[0]=sum/m;

for(suni=0,i=0;i<m;i++){

sum+=x[i];

}

alf[l]=sum/m;

n[l]=l；

Pfl][l]=l；

p[l][O]=-alfTl];

for(sum=0,i=0;i<iTi;i++){

sum+=y[i]*(x[i]-alf[1J);

}

for(sum1=0,i=0;i<ni;i++){

suml+=(x[i]-alf[l])*(x[i]-alf[l]);

)

a[l]=sum/suml;

addpoly(pf01,n[0],a[01,p[l],n[l],a[l],result,&nr);

k=l;

do{

k++;

for(sum=0,i=0;i<m;i++){

temp=Po