數(shù)值分析實(shí)驗(yàn)程序

數(shù)值分析實(shí)驗(yàn)程序...[推薦]

1）Guass列主元消元法解線性方程組：

//高斯列主元消去法求解非線性方程組源程序清單

#defineNUMBER20

#include"math.h"

floatA[NUMBER][NUMBER+1],ark;

intflag,n;

exchange(intr,intk)

{inti;

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

A[0][i]=A[r][i];

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

A[r][i]=A[k][i];

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

A[k][i]=A[0][i];

}

floatmax(intk)

{inti;

floattemp=0;

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

if(fabs(A[i][k])>temp)

{temp=fabs(A[i][k]);

flag=i;}

returntemp;

}

main()

{floatx[NUMBER],b[NUMBER];

intr,k,i,j;

me:printf("\n\n\n用Gauss列主元消元法解線性方程組");

printf("\ninputn=");

scanf("%d",&n);

printf("\n\n現(xiàn)在輸入系數(shù)矩陣A和向量b:\n\n");

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

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

{if(j==n+1)printf("請(qǐng)輸入b%d:",i);

elseprintf("請(qǐng)輸入a%d%d:",i,j);

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

}

k=1;

label:

/****************************************************/

ark=max(k);

/*if(fabs(A[1][k])>=fabs(A[2][k]))

{ark=fabs(A[1][k]);flag=1;}

else{ark=fabs(A[2][k]);flag=2;}

if(fabs(A[3][k])>ark)

{ark=fabs(A[3][k]);flag=3;}*/

/**************************************************/

if(ark==0)

{printf("\n不是合法的方程組!");exit(0);}

elseif(flag!=k)

exchange(flag,k);

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

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

A[i][j]=A[i][j]-A[i][k]*A[k][j]/A[k][k];

if(k<n-1)

{k++;gotolabel;}

else

{x[n]=A[n][n+1]/A[n][n];

for(k=n-1;k>=1;k--)

{floatme=0;

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

{me=me+A[k][j]*x[j];}

x[k]=(A[k][n+1]-me)/A[k][k];

}

}

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

printf("\nx%d=%f",i,x[i]);

gotome;

}

///////////////////////////////////////////////////////////////////////

2）//插值源程序

#definef(x)1/(1+25*x*x)

#definef1(x)-50*x/((1+25*x*x)*(1+25*x*x))

#include"math.h"

floatxi[11],M[11];

intN=10;

floath=0.2;

floatL(inti,floatx)

{floatll=1.0;intj;

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

{if(i==j)continue;

elsell=ll*(x-xi[j]);

}

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

{if(i==j)continue;

elsell=ll/(xi[i]-xi[j]);

}

returnll;

}

floatS(intN,floatx)

{floatss=0;intj;

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

ss=ss+f(xi[j])*L(j,x);

returnss;

}

floatS1(floatx0,floatx1,floatx)

{floatss;

ss=(x-x1)/(x0-x1)*f(x0)+(x-x0)/(x1-x0)*f(x1);

returnss;

}

floatchashan(floatx0,floatx1,floatx2)

{return((f(x2)-f(x0))/(x2-x0)-(f(x1)-f(x0))/(x1-x0))/(x2-x1);

}

floatS3(floatxj,floatxjp1,floatMj,floatMjp1,floatx)

{floatss1,ss2,ss3,ss4,ss;

ss1=(xjp1-x)*(xjp1-x)*(xjp1-x)/(6*h)*Mj;

ss2=(x-xj)*(x-xj)*(x-xj)/(6*h)*Mjp1;

ss3=(x-xj)/h*(f(xjp1)-h*h/6*Mjp1);

ss4=(xjp1-x)/h*(f(xj)-h*h/6*Mj);

ss=ss1+ss2+ss3+ss4;

returnss;

}

main()

{inti,j;intcontrol;floathelp1,help2,helpe,helpb;

floatd[11];intcchar=0;

floatx0,xo1[10],xo2[10];

x0=0;

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

xi[i]=-1.0+i*2.0/N;

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

{xo1[j]=0.06+0.1*j;

xo2[j]=-0.06-0.1*j;

}

label:

printf("\n\nRunge現(xiàn)象的克服");

printf("\n\n1.多項(xiàng)式插值");

printf("2.分段線性插值");

printf("3.三次樣條函數(shù)插值");

printf("4.退出");

printf("\n請(qǐng)選擇:");

scanf("%d",&control);

switch(control)

{case1:

//下面是多項(xiàng)式插值

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

printf("\nf(%5.2f)=%9.6f",xo1[i],f(xo1[i]));

printf("SN(%5.2f)=%9.6f",xo1[i],S(N,xo1[i]));

printf("\nf(%5.2f)=%9.6f",xo2[i],f(xo2[i]));

printf("SN(%5.2f)=%9.6f",xo2[i],S(N,xo2[i]));

help1=f(xo1[i]);help2=S(N,xo1[i]);

helpe=help1-help2;

printf("f-SN=e=%9.6f",helpe);

helpb=helpe/help1*100.0;

printf("相對(duì)誤差=%9.6f%%",helpb);

}

break;

case2:

//下面是分段線性插值

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

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

{if(xo1[j]>xi[i]&&xo1[j]<xi[i+1])

{help1=f(xo1[j]);

help2=S1(xi[i],xi[i+1],xo1[j]);

helpe=help1-help2;

helpb=helpe/help1*100.0;

printf("\nf(%5.2f)=%9.6f",xo1[j],help1);

printf("S1(%5.2f)=%9.6f",xo1[j],help2);

printf("f-S1=e=%9.6f",helpe);

printf("相對(duì)誤差=%9.6f%%",helpb);

}

if(xo2[j]>xi[i]&&xo2[j]<xi[i+1])

{help1=f(xo2[j]);

help2=S1(xi[i],xi[i+1],xo2[j]);

printf("\nf(%5.2f)=%9.6f",xo2[j],help1);

printf("S1(%5.2f)=%9.6f",xo2[j],help2);

}

}

break;

case3:

//以下是三次樣條函數(shù)插值

d[0]=6/h*((f(xi[1])-f(xi[0]))/h-f1(xi[0]));

printf("\nd0=%f",d[0]);

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

{d[j]=6*chashan(xi[j-1],xi[j],xi[j+1]);

printf("\nd%d=%f",j,d[j]);

}

d[N]=6/h*(f1(xi[N])-(f(xi[N])-f(xi[N-1]))/h);

printf("\nd10=%f",d[N]);

printf("\n\n請(qǐng)用Gauss列主元法解十一階三對(duì)角方程求出M0至M10

!");

printf("\n然后把M0至M10輸入!準(zhǔn)備好后按1:");

scanf("%d",&cchar);

if(cchar==1)

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

{printf("\n請(qǐng)輸入M%d:",j);

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

}}

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

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

{if(xo1[j]>xi[i]&&xo1[j]<xi[i+1])

{help1=f(xo1[j]);

help2=S3(xi[i],xi[i+1],M[i],M[i+1],xo1[j]);

helpe=help1-help2;

helpb=helpe/help1*100.0;

printf("\nf(%5.2f)=%9.6f",xo1[j],help1);

printf("S3(%5.2f)=%9.6f",xo1[j],help2);

printf("f-S1=e=%9.6f",helpe);

printf("相對(duì)誤差=%9.6f%%",helpb);

}

if(xo2[j]>xi[i]&&xo2[j]<xi[i+1])

{help1=f(xo2[j]);

help2=S3(xi[i],xi[i+1],M[i],M[i+1],xo2[j]);

printf("\nf(%5.2f)=%9.6f",xo2[j],help1);

printf("S3(%5.2f)=%9.6f",xo2[j],help2);

}

}

break;

case4:exit(0);

break;

}

gotolabel;

}



//////////////////////////////////////////////////////////////////////////////

/

3）//二分法求解非線性方程

#include"math.h"

#definef(x)x*x*x-x-1

main()

{floata,b,c,e;

intk=0;

charaa;

printf("\n\n\n二分法求非線性方程的解");

printf("\n\n為使程序簡(jiǎn)單，f(x)在源程序中定義!");

printf("\n已知f(x)在[a,b]上連續(xù)");

printf("\n請(qǐng)輸入a:");

scanf("%f",&a);

scanf("%c",&aa);

printf("請(qǐng)輸入b:");

scanf("%f",&b);

scanf("%c",&aa);

printf("請(qǐng)輸入允許誤差e:");

scanf("%f",&e);

scanf("%c",&aa);

label:

c=(a+b)/2;

printf("\n\na%d=%f",k,a);

printf("\nb%d=%f",k,b);

printf("\nx%d=%f",k,c);

k++;

if((f(a))*(f(c))>0)a=c;

elseb=c;

if(b-a<=e)

{printf("\n最終結(jié)果是:");

printf("\nx%d=%f",k,(a+b)/2);

exit(0);

}

elsegotolabel;}

//////////////////////////////////////////////////////////////////////////////

/

4）//牛頓迭代法求非線性方程的解源程序清單

#include"math.h"

#definef(x)x*x*x-x-1

#definef1(x)3*x*x-1

main()

{floata,x0,x1,e1,e2;

intN,k;charaa;char*str1,*str2;

printf("\n\n\n用牛頓迭代法求非線性方程的解");

printf("\n\n\n為使程序簡(jiǎn)單，函數(shù)f(x)及f'(x)在源程序中定義!");

printf("\n\n請(qǐng)輸入初值x0:");

scanf("%f",&x0);

/*scanf("%c",&aa);*/

printf("輸入誤差e1:");

scanf("%f",&e1);

scanf("%c",&aa);

printf("輸入誤差e2:");

scanf("%f",&e2);

scanf("%c",&aa);

printf("輸入最大迭代次數(shù)N:");

scanf("%d",&N);

scanf("%c",&aa);

k=1;

label:

if(f(x0)==0){printf("\n%f為所求得根之一!",x0);

exit(0);}

elseif(fabs(f(x0))<e1){printf("\n奇異標(biāo)志!建議重取初值或取更小的e1,e2!")

;

exit(0);}

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

printf("\n第%d次迭代:",k);

printf("\nx%d=%f",k,x1);

if(fabs(x1-x0)<e2){printf("\n%f為所求近似解!",x1);

exit(0);}

elseif(k==N){printf("\n迭代失敗!建議取更大的N值或另選初值!");exit

(0);}

else{x0=x1;

k++;

}

gotolabel;

}//////////////////////////////////////////////////////////////////////////////

/

5）//Romberg積分法

#include"math.h"

floatT[20][20];

floatf(floatx)

{floatss;

if(x==0)x=1e-5;//防止除數(shù)為零而溢出！

ss=sin(x*x);//為被積函數(shù)

returnss;

}

main()

{floata,b,e;floatsum=0.0;floattemp,help;

inti,j,k,m,p,q;

printf("\nRomberg