数值分析计算方法程序汇总(2)

y[i+1]=exp(x[i+1]);

p[i+1]=y[i]+(y[i+1]-y[i])/(x[i+1]-x[i])*(xl-x[i]);

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

}

结果：输入i=0时，p[i+1]=3.214028

（十五）解线性方程组

5x1 x2 x3 1

例：解方程组x1 4x2 x3 2

x1 x2 6x3 3

程序：

#include "stdio.h"

void main()

{double a[3][3]={{5.0,1.0,1.0},{1.0,4.0,-1.0},{1.0,-1.0,6.0}},b[3]={1.0,2.0,3.0},x[3]; int i,j,k,n;

n=3;

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

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

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

{a[i][j]=a[i][j]-a[k][j]*a[i][k]/a[k][k];

b[i]=b[i]-b[k]*a[i][k]/a[k][k];

}

}

}

x[n-1]=b[n-1]/a[n-1][n-1];

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

{x[k]=(b[k]-a[k][k]*x[k+1]-a[k][n]*x[n])/a[k][k];

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

}

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

}

结果：-0.414921

0.414921

0.572816

（十六）列主元消去法

x1 x2 x3 3

例：解方程组x1 x2 x3 1

2x1 x2 x3 5

程序：

#include<stdio.h>

#include<math.h>

#define N 100

float a[N][N+1];

void main( )

{ int i,j,k,n;

float t,s=0;

printf("输入矩阵阶数:"); scanf("%d",&n);

printf("\n"); printf("输入增广矩阵:\n");

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

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

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

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

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

if( fabs(a[i][k]) > fabs(a[k][k]) )

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

{ t=a[k][j]; a[k][j]=a[i][j]; a[i][j]=t;}

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

{ a[i][k]=a[i][k] / a[k][k];

for(j=k+1;j<n+1;j++) a[i][j]=a[i][j]-a[i][k] * a[k][j];

}

}

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

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

{s=0;

for(j=k+1;j<n;j++) s+=a[k][j]*a[j][n]; a[k][n]=( a[k][n]-s ) / a[k][k]; }

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

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

}

结果：x[1]=2.000000

x[2]=1.000000

x[3]=-0.000000

（十七）追赶法

4x1 x2 1

例：x1 4x2 x3 2

x2 4x3 x4 1

x3 4x4 0

程序：

#include "stdio.h"

void main()

{double x[3],a[4]={0,1,1,1},b[4]={4,4,4,4},c[4]={1,1,1,0},d[4]={1,2,1,0},e[4],f[4]; int k,n=3;

e[0]=d[0]/b[0];f[0]=c[0]/b[0];

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

{e[k]=(d[k]-a[k]*e[k-1])/(a[k]*f[k-1]+b[k]);

f[k]=c[k]/(a[k]*f[k-1]+b[k]);

}

x[n]=(d[n]-a[n]*e[n-1])/(b[n]-a[n]*f[n-1]);

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

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

{x[k]=e[k]-f[k]*x[k+1];

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

}

}

结果：-0.036900

0.147601

0.377035

0.155741

（十八）高斯-赛德尔迭代法

8x1 3x2 2x3 20

例：4x1 11x2 x3 33

6x1 3x2 12x3 36

程序：

#include "stdio.h"

#include <math.h>

#define N 100

void main()

{double a[N][N]={{8,-3,2},{4,11,-1},{6,3,12}},b[N]={20,33,36},x[N]={0,0,0},xl,s; int i,j,k,n=3;

for(k=1;k<=10;k++)

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

{s=0.0;

xl=x[i];

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

{if(j==i)

continue;

else s=s+a[i][j]*x[j];

}

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

}

}

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

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

}

结果：

x[1]=3.000000

x[2]=2.000000

x[3]=1.000000

（十九)一步欧拉公式

yk 1 yk hf(xk,yk)

例：已知y' x y

y(0) 1求x=1时的值。

程序：

#include "stdio.h"

void main()

{

double x,y,h;

int i,n=10;

h=0.1;x=0;y=1;

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

{

y=y+h*(x+y);

x=x+h;

}

printf("%f\n",y);

}

结果：3.187485

（二十）隐式欧拉公式：

yk 1 yk hf(xk 1,yk 1)

例：已知y' x y

y(0) 1求x=1时的值。

程序：

#include "stdio.h"

#include <math.h>

void main()

{

double x[10],y[10],yl,h,e;

int i,j,n=10;

h=0.1;e=0.01;

x[0]=0.0;y[0]=1.0;

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

{

y[i+1]=y[i]+h*(x[i]+y[i]);

x[i+1]=x[i]+h;

for(j=0;;j++)

{

y[j+1]=y[j]+h*(x[j+1]+y[j+1]);

yl=y[j+1];

if(fabs(yl-y[i+1])<e)

break;

else y[j+1]=yl;

}

}

printf("%f\n",yl);

}

结果：3.728783

（二十一）四阶龙格-库塔公式

例：已知y'=x+y，y(0)=1；求y(0.2)=？

程序：

#include "stdio.h"

void main()

{double x,y,h,k1,k2,k3,k4;

int i,n;

h=0.1;n=2;x=0.0;y=1.0;

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

{

k1=x+y;

k2=x+y+h/2+h*k1/2;

k3=x+y+h/2+h*k2/2;

k4=x+y+h+h*k3;

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

x=x+h;

}

printf("y(0.2)=%f\n",y);

}

结果：y(0.2)=1.242805

补充例题：

（二十二）

例：已知y''=x+y+y'，用一步欧拉公式求n=10时，y=？，y'=？ 程序：

#include "stdio.h"

void main()

{

double x,y,z,h;

int i,n;

h=0.1;n=10;

x=0;y=0;z=1;

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

{

y=y+h*z;

z=z+h*(x+y+z);

x=x+h;

}

printf("y=%f,y'=%f\n",y,z);

}

结果：y=2.001568,y'=4.314725

（二十三）

例：已知y''=x+y，y(0)=1，y(1)=0，求n=3时，y(1)=？y(2)=？ 程序：

#include "stdio.h"

void main()

{

double x[3],y[3],d[3],e[3],f[3],b,h;

int i,k,n=3;

h=1.0/n;b=-2.0-h*h;

x[1]=h;d[1]=x[1]*h*h-1.0;

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

{

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

d[i]=h*h*x[i];

}

e[1]=d[1]/b;f[1]=-1/b;

for(k=2;k<=n-1;k++)

{

e[k]=(d[k]-e[k-1])/(b+f[k-1]);

f[k]=-1.0/(b+f[k-1]);

}

y[n-1]=e[n-1];

printf("y[%d]=%f",n-1,y[n-1]);

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

{

y[k]=e[k]+f[k]*y[k+1]; printf("y[%d]=%f",k,y …… 此处隐藏：1915字，全部文档内容请下载后查看。喜欢就下载吧 ……