实验一：拉格朗日插值多项式
命名（源程序.cpp及工作区.dsw）：lagrange
问题：
4
//Lagrange.cpp
#include
#include
#define N 4
int checkvalid(double x[], int n);
void printLag (double x[], double y[], double varx, int n);
double Lagrange(double x[], double y[], double varx, int n);
void main ()
{
double x[N+1] = {0.4, 0.55, 0.8, 0.9, 1};
double y[N+1] = {0.41075, 0.57815, 0.88811, 1.02652, 1.17520};
double varx = 0.5;
if (checkvalid(x, N) == 1)
{
printf("\n\n插值结果: P(%f)=%f\n", varx, Lagrange(x, y, varx, N));
}
else
{
printf("结点必须互异");
}
getch();
}
int checkvalid (double x[], int n)
{
int i,j;
for (i = 0; i < n; i++)
{
for (j = i + 1; j < n+1; j++)
{
if (x[i] == x[j])//若出现两个相同的结点，返回-1
{
return -1;
}
}
}
return 1;
}
double Lagrange (double x[], double y[], double varx, int n) {
double fenmu;
double fenzi;
double result = 0;
int i,j;
printf("Ln(x) =\n");
for (i = 0; i < n+1; i++)
{
fenmu = 1;
for (j = 0; j < n+1; j++)
{
if (i != j)
{
fenmu = fenmu * (x[i] - x[j]);
}
}
printf("\t%f", y[i] / fenmu);
fenzi = 1;
for (j = 0; j < n+1; j++)
{
if (i != j)
{
printf("*(x-%f)", x[j]);
fenzi = fenzi * (varx - x[j]);
}
}
if (i != n)
{
printf("+\n");
}
result += y[i] / fenmu * fenzi;
}
return result;
}运行及结果显示：
实验二：牛顿插值多项式
命名（源程序.cpp及工作区.dsw）：newton_cz
问题：
4
//newton_cz.cpp
#include
#include
#include
#define N 4
int checkvaild(double x[],int n)
{
int i,j;
for(i=0;i { for(j=i+1;j if(x[i]==x[j]) return -1; } return 1; } void chashang(double x[],double y[],double f[][N+1]) { int i,j,t=0; for(i=0;i { f[i][0]=y[i];//f[0][0]=y[0],f[1][0]=y[1];f[2][0]=y[2];f[3][0]=y[3];f[4][0]=y[4] } for(j=1;j { for(i=0;i f[i][j]=(f[i+1][j-1]-f[i][j-1])/(x[t+i+1]-x[i]);//一阶为f[0][1]~f[3][1];二阶为f[0][2]~f[2][2] //三阶为f[0][3]~f[1][3];四阶为f[0][4] t++; } } double compvalue(double t[][N+1],double x[],double varx) { int i,j,r=0; double sum=t[0][0],m[N+1]={sum,1,1,1,1}; printf("i\tXi\t F(Xi)\t 1阶\t 2阶\t\t3阶\t 4阶\n"); printf("--------------------------------------------------------------------------------"); for(i=0;i { r=i; printf("x%d\t%f ",i,x[i]); for(j=0;j<=i;j++) { printf("%f ",t[r][j]); r--; } printf("\n"); } printf("--------------------------------------------------------------------------------");/**/ printf("N(x) =\n"); printf(" %f\n",t[0][0]); for(i=1;i { for(j=0;j m[i]=m[i]*(varx-x[j]); m[i]=t[0][i]*m[i]; sum=sum+m[i]; } for(i=1;i { printf(" +%f*",t[0][i]); for(j=0;j printf("(x-%f)",x[j]); printf("\n"); } return sum; } void main() { double varx,f[N+1][N+1]; double x[N+1]={0.4,0.55,0.8,0.9,1}; double y[N+1]={0.41075,0.57815,0.88811,1.02652,1.17520}; checkvaild(x,N); chashang(x,y,f); varx=0.5; if(checkvaild(x,N)==1) { chashang(x,y,f); printf("\n\n牛顿插值结果: P(%f)=%f\n",varx,compvalue(f,x,varx)); } else printf("输入的插值节点的x值必须互异!"); getch(); } 运行及结果显示： 实验三：自适应梯形公式 命名（源程序.cpp 及工作区.dsw ）：autotrap ][n T 问题：计算? +=1 02 14dx x π的近似值，使得误差（EPS ）不超过6 10- //autotrap.cpp #include return 4/(1+x*x); } double AutoTrap(double(*f)(double),double a,double b,double eps) { double t2,t1,sum=0.0; int i,k; t1=(b-a)*(f(a)+f(b))/2; printf("T(%4d)=%f\n",1,t1); for(k=1;k<11;k++) { for(i=1;i<=pow(2,k-1);i++) sum+=f(a+(2*i-1)*(b-a)/pow(2,k)); t2=t1/2+(b-a)*sum/pow(2,k); printf("T(%4d)=%f\n",(int)pow(2,k),t2); if(fabs(t2-t1) return sum; } void main() { double s; s=AutoTrap(f,0.0,1.0,EPS); getch(); } 运行及结果显示： 实验四：龙贝格算法 命名（源程序.cpp 及工作区.dsw ）：romberg 问题：求? +=1 02 14dx x π的近似值，要求误差不超过6 10- //romberg.cpp #include double f(double x){return 4/(1+x*x);} double Romberg(double a,double b,double (*f)(double),double eps) { double T[N][N],p,h; int k=1,i,j,m=0; T[0][0]=(b-a)/2*(f(a)+f(b)); do{ p=0; h=(b-a)/pow(2,k-1); for(i=1;i<=pow(2,k-1);i++) p=p+f(a+(2*i-1)*h/2); T[0][k]=T[0][k-1]/2+p*h/2; for(i=1;i<=k;i++) {j=k-i; T[i][j]=(pow(4,i)*T[i-1][j+1]-T[i-1][j])/(pow(4,i)-1); } k++; p=fabs(T[k-1][0]-T[k-2][0]); }while(p>=EPS); k--; while(m<=k) {for(i=0;i<=m;i++) printf("T(%d,%d)=%f ",i,m-i,T[i][m-i]); m++; printf("\n"); } return T[k][0]; } void main() { printf("\n\n 积分结果 = %f",Romberg(0,1,f,EPS)); getch(); } 运行及结果显示： 实验五：牛顿切线法 问题：求方程01)(3=--=x x x f 在5.1=x ，6.0=x 附近的根（精度3102 1 -?= ） //newton_qxf.cpp #include double f(double x) { return x*x*x-x-1; } double f1(double x) { return 3*x*x-1; } int newton(double (*f)(double), double (*f1)(double), double &x0, double eps) { double xk, xk1; int count = 0; printf("k\txk\n"); printf("-----------------------\n"); xk = x0; printf("%d\t%f\n", count, xk); do { if((*f1)(xk)==0) return 2; xk1 = xk - (*f)(xk) / (*f1)(xk); if (fabs(xk - xk1) < eps) { count++; xk = xk1; printf("%d\t%f\n", count, xk); x0 = xk1; return 1; } count++; xk = xk1; printf("%d\t%f\n", count, xk); }while(count < MaxK); return 0; } void main() { for(int i=0;i<2;i++) { double x=0.6; if(i==1) x=1.5; printf("x0初值为%f:\n",x); if (newton(f, f1, x, EPS) == 1) { printf("-----------------------\n"); printf("the root is x=%f\n\n\n", x); } else { printf("the method is fail!"); } } getch(); } 实验六：牛顿下山法 命名（源程序.cpp 及工作区.dsw ）：newton_downhill 问题：求方程01)(3=--=x x x f 在5.1=x ，6.0=x 附近的根（精度3102 1 -?=） //newton_downhill.cpp #include #define Et 1e-3//下山因子下界 #define E1 1e-3//根的误差限 #define E2 1e-3//残量精度 double f(double x) { return x * x * x - x - 1; } double f1(double x) { return 3 * x * x - 1; } void errormess(int b) { char *mess; switch(b) { case -1: mess = "f(x)的导数为0!"; break; case -2: mess = "下山因子已越界，下山处理失败"; break; default: mess = "其他类型错误!"; } printf("the method has faild! because %s", mess); } int Newton(double (*f)(double), double (*f1)(double), double &x0) { int k = 0; double t; double xk, xk1; double fxk, fxk_temp; printf("k t xk f(xk)\n"); printf("----------------------------------------------------------\n"); printf("%-20d", k); xk = x0; printf("%-15f", x0); fxk = (*f)(xk); printf("%-20f", fxk); printf("\n"); for (k = 1; ; k++) { t = 1; while(1) { printf("%-10d", k); printf("%-10f", t); if((*f1)(xk) != 0) { xk1 = xk - t * (*f)(xk) / (*f1)(xk); } else { return -1; } printf("%-15f", xk1); fxk_temp = (*f)(xk1); printf("%-20f", fxk_temp); if(fabs(fxk_temp) >fabs(fxk)) { t = t / 2; printf("\n"); if (t < Et) { return -2; } } else { printf("下山成功\n"); break; } } if (fabs(xk-xk1) { x0 = xk1; return 1; } xk = xk1; } } void main() { int b; double x0; x0 = 0.6; b = Newton(f, f1, x0); if (b == 1) printf("\nthe root x=%f\n", x0); else errormess(b); getch(); } 运行及结果显示： 实验七：埃特金加速算法 命名（源程序.cpp 及工作区.dsw ）：aitken 问题：求方程01)(3=--=x x x f 在5.1=x ，6.0=x 附近的根（精度3102 1 -?= ） //aitken.cpp #include double g(double x) { return x * x * x - 1; } int aitken (double (*g)(double), double &x, double eps) { int k; double xk = x, yk, zk, xk1; printf("k xk yk zk xk+1\n"); printf("-------------------------------------------------------------------\n"); for (k = 0;k xk1 = xk - (yk - xk) * (yk - xk) / (zk - 2 * yk + xk); printf("%-10d%-15f%-15f%-15f%-15f\n", k, xk, yk, zk, xk1); if (fabs(yk-xk)<=eps) { x = xk1; return k + 1; } xk = xk1; } return -1; } void main () { double x = 1.5; int k; k = aitken(g, x, EPS); if (k == -1) printf("迭代次数越界!\n"); else printf("\n 经k=%d 次迭代，所得方程根为:x=%f\n", k, x); getch(); } 运行及结果显示： 实验八：正割法 问题：求方程01)(3=--=x x x f 在5.1=x 附近的根（精度0.5e-8） //ZhengGe.cpp #include #define MaxK 1000 #define EPS 0.5e-8 double f(double x) { return x*x*x-x-1; } int ZhengGe(double (*f)(double), double x0, double &x1, double eps) { printf("k xk f(xk)\n"); printf("---------------------------------------------\n"); int k; double xk0, xk, xk1; xk0 = x0; printf("%-10d%-15f%-15f\n", 0, x0, (*f)(x0)); xk = x1; for (k=1;k { if((*f)(xk)-(*f)(xk0)==0) return -1; xk1 = xk - (*f)(xk) / ( (*f)(xk) - (*f)(xk0) ) * (xk - xk0); printf("%-10d%-15f%-15f\n", k, xk, (*f)(xk)); if (fabs(xk1 - xk)<=eps) { printf("%-10d%-15f%-15f\n\n", k + 1, xk1, (*f)(xk1)); printf("---------------------------------------------\n\n"); x1 = xk1; return 1; } xk0 = xk; xk = xk1; } return -2; } void main () { double x0 = 1, x1 = 2; if (ZhengGe(f, x0, x1, EPS) == 1) { printf("the root is x = %f\n", x1); } else { printf("the method is fail!"); } getch(); } 实验九：高斯列选主元算法 命名（源程序.cpp 及工作区.dsw ）：colpivot 问题：求解方程组并计算系数矩阵行列式值 ??? ??=-+=+-=-+2 24053 2321 321321x x x x x x x x x //colpivot.cpp #include static double aa[N][N]={{1,2,-1},{1,-1,5},{4,1,-2}}; static double bb[N+1]={3,0,2}; void main() { int i,j; double a[N+1][N+1],b[N+1],x[N+1],det; double gaussl(double a[][N+1],double b[],double x[]); for(i=1;i<=N;i++) { for(j=1;j<=N;j++) a[i][j]=aa[i-1][j-1]; b[i]=bb[i-1]; } det=gaussl(a,b,x); if(det!=0) { printf("\n 方程组的解为:"); for(i=1;i<=N;i++) printf(" x[%d]=%f",i,x[i]); printf("\n\n 系数矩阵的行列式值＝%f",det); } else printf("\n\n 系数矩阵奇异阵，高斯方法失败 !:"); getch(); } double gaussl(double a[][N+1],double b[],double x[]) //a传入增广矩阵(有效元素为a[1,1]...a[n,n+1]),x负责传入迭代初值并返回解向量//(有效值为x[1]...x[n]);返回值：系数矩阵行列式值detA { double det=1.0,F,m,temp; int i,j,k,r; void disp(double a[][N+1],double b[]); printf("消元前增广矩阵:\n"); disp(a,b); for(k=1;k { temp=a[k][k]; r=k; for(i=k+1;i<=N;i++) { if(fabs(temp) { temp=a[i][k]; r=i; }//按列选主元，即确定ik } if(a[r][k]==0) return 0;//如果aik,k=0,则A为奇异矩阵，停止计算 if(r!=k) { for(j=k;j<=N;j++) { a[k][j]+=a[r][j]; a[r][j]=a[k][j]-a[r][j]; a[k][j]-=a[r][j]; } b[k]+=b[r]; b[r]=b[k]-b[r]; b[k]-=b[r]; det=-det;//如果ik!=k,则交换[A,b]第ik行与第k行元素 printf("交换第%d, %d行:\n",k,r); disp(a,b); } for(i=k+1;i<=N;i++) { m=a[i][k]/a[k][k]; for(j=1;j<=k;j++) a[i][j]=0.0; for(j=k+1;j<=N;j++) a[i][j]-=m*a[k][j]; b[i]-=m*b[k]; } printf("第%d次消元:\n",k); disp(a,b); det*=a[k][k]; } //大FOR循环结束 x[N]=b[N]/a[N][N]; for(i=N-1;i>0;i--) { F=0.0; for(j=i+1;j<=N;j++) F+=a[i][j]*x[j]; x[i]=(b[i]-F)/a[i][i]; } det*=a[N][N]; return det; } void disp(double a[][N+1],double b[]) //显示选主元及消元运算中间增广矩阵结果 {int i,j; for(i=1;i<=N;i++) {for(j=1;j<=N;j++) printf("%10f\t",a[i][j]); printf("%10f\n",b[i]); } } 运行及结果显示： 实验十：高斯全主元消去算法 命名（源程序.cpp 及工作区.dsw ）：fullpivot 问题：利用全主元消去法求解方程组?? ?? ??????-=????????????????????----9555.04525.015.0201 .0321x x x //fullpivot.cpp #include static double aa[N][N]={{0.01,2,-0.5},{-1,-0.5,2},{5,-4,0.5}}; static double bb[N]={-5,5,9}; void main() { int i,j; double a[N+1][N+1],b[N+1],x[N+1],det; int t[N+1];//引入列交换保存x向量的各分量的位置顺序 double gaussl(double a[][N+1],double b[],double x[],int t[]); for(i=1;i<=N;i++) { for(j=1;j<=N;j++) a[i][j]=aa[i-1][j-1]; b[i]=bb[i-1]; } for(i=1;i<=N;i++) t[i]=i; det=gaussl(a,b,x,t); if(det!=0) { printf("\n方程组的解为:"); for(i=N;i>=1;i--) { printf(" x[%d]=%f",t[i],x[i]); if(i>1)printf(" -->"); } printf("\n\n系数矩阵的行列式值＝%f",det); } else printf("\n\n系数矩阵奇异阵，高斯方法失败!:"); getch(); } double gaussl(double a[][N+1],double b[],double x[],int t[]) //a传入增广矩阵(有效元素为a[1,1]...a[n,n+1]),x负责传入迭代初值并返回解向量//(有效值为x[1]...x[n]);返回值：系数矩阵行列式值detA { double det=1.0,F,m,temp; int i,j,k,r,s; void disp(double a[][N+1],double b[],int x[]); // printf("消元前增广矩阵:\n"); //disp(a,b,t); for(k=1;k { temp=a[k][k]; r=k; s=k; for(i=k;i<=N;i++) for(j=k;j<=N;j++) if(fabs(temp) { temp=a[i][j]; r=i;s=j; }//选主元，选取ik,jk if(a[r][s]==0) return 0; if(r!=k) { for(j=k;j<=N;j++) { a[k][j]+=a[r][j]; a[r][j]=a[k][j]-a[r][j]; a[k][j]-=a[r][j]; } b[k]+=b[r]; b[r]=b[k]-b[r]; b[k]-=b[r]; det=-det; printf("交换第%d, %d行:\n",k,r); disp(a,b,t); } if(s!=k) { for(i=0;i<=N;i++) { a[i][k]+=a[i][s]; a[i][s]=a[i][k]-a[i][s]; a[i][k]-=a[i][s]; } t[k]+=t[s]; t[s]=t[k]=t[s]; t[k]-=t[s]; det=-det; printf("交换第%d, %d列:\n",k,s); disp(a,b,t); } for(i=k+1;i<=N;i++) { m=a[i][k]/a[k][k]; for(j=1;j<=k;j++) a[i][j]=0.0; for(j=k+1;j<=N;j++) a[i][j]-=m*a[k][j]; b[i]-=m*b[k]; } //消元计算printf("第%d次消元:\n",k); disp(a,b,t); det*=a[k][k]; } //大FOR循环结束 x[N]=b[N]/a[N][N]; for(i=N-1;i>0;i--) { F=0.0; for(j=i+1;j<=N;j++) F+=a[i][j]*x[j]; x[i]=(b[i]-F)/a[i][i]; }//回代计算 det*=a[N][N]; return det; } void disp(double a[][N+1],double b[],int x[]) //显示选主元及消元运算中间增广矩阵结果 { int i,j; for(i=1;i<=N;i++) { for(j=1;j<=N;j++) printf("%f\t",a[i][j]); printf("| x%d = %f\n",x[i],b[i]); } } 运行及结果显示： 实验十一：LU 分解算法 命名（源程序.cpp 及工作区.dsw ）：LU 问题：利用LU 法求解方程组 ?? ?? ? ?????=????????????????????472111312613321x x x //LU.cpp #include static aa[N][N]={{3,1,6},{2,1,3},{1,1,1}}; static bb[N]={2,7,4}; void main() { int i,j; double a[N+1][N+1],b[N+1]; void solve(double a[][N+1],double b[N+1]); int LU(double a[][N+1]); for(i=1;i<=N;i++) { for(j=1;j<=N;j++) a[i][j]=aa[i-1][j-1]; b[i]=bb[i-1]; } if(LU(a)==1) 第2章矩阵与线性代数方程组 一般的线性代数方程组，A非奇异可根据Cramer法则求解方程唯一解但是它的计算量很大。 高斯消元法的算法时间复杂度是O(n3)，可以解一系列的线性方程；所占数据空间符合原地工作的原则。但是算法对数值计算不稳定（当分母为0或很小时）。可以用在计算机中来解决数千条等式及未知数。不过，如果有过百万条等式时，这个算法会十分费时。 解决高斯法中的不稳定性，在每次归一化前增加选主元（列选主元、全选主元）过程。但是列选主元法仍不稳定，不适求解大规模线性代数方程组。全选主元的高斯消去法，则在复杂度降低的同时能够避免舍入误差，保证数值稳定性。 高斯-约当消去法算法产生出来的矩阵是一个简化行梯阵式，而不是高斯消元法中的行梯阵式。相比起高斯消元法，此算法的效率比较低，却可把方程组的解用矩阵一次过表示出来。线性代数方程组的迭代解法 简单迭代法：迭代格式发散但迭代值序列不一定发散，但收敛格式收敛，迭代值序列收敛于方程组的准确解与选取迭代初值无关。 雅可比迭代法: 计算公式简单，且计算过程中原始矩阵A始终不变，比较容易并行计算。但是收敛速度较慢，而且占据的存储空间较大，所以工程中一般不直接用雅克比迭代法，而用其改进方法。 高斯-赛德尔迭代法：较上面的迭代复杂，但是矩阵的条件相对宽松。 松弛法：需要根据经验去调整，收敛速度依赖松弛参数的选择，收敛条件的要求更宽松。共轭梯度法：是介于最速下降法与牛顿法之间的一个方法，它仅需利用一阶导数信息，但克服了最速下降法收敛慢的缺点，又避免了牛顿法需要存储和计算Hesse矩阵并求逆的缺点，在各种优化算法中，共轭梯度法是非常重要的一种。其优点是所需存储量小，具有步收敛性，稳定性高，而且不需要任何外来参数。共轭梯度法不仅是解决大型线性方程组最有用的方法之一，也是解大型非线性最优化最有效的算法之一。 第3章矩阵特征值 乘幂法计算绝对值最大的特征值：其收敛速度受限于最大与次大特征值比值绝对值的大小，实际应用中采用加速技术。 求对称特征值的雅克比方法96：每进行一次选装变换钱都需要在飞对角线的元素中选取绝对值最大的元素，很费时间，雅克比过关法对此做了改进。 QR方法求一般实矩阵的全部特征值98下100下：重复多次进行QR分解费时，计算工作量很大。一般先进行相似变换然后进行QR分解。但是这样仍然收敛速度慢，一般是线性收敛。实际应用中使用双重步QR变换将带原点的QR算法中相邻两步合并一步，加速收敛避免复数运算。 第4章非线性方程与方程组 二分法：每次运算后，区间长度减少一半，是线形收敛。优点是简单，但是不能计算复根和重根。 简单迭代法：直接的方法从原方程中隐含的求出x，从而确定迭代函数 (x)，这种迭代法收敛速度较慢，迭代次数多。 埃特金迭代法113中：对简单迭代进行改进，使在其不满足收敛条件下迭代过程也收敛，在其收敛时加快收敛速度，减少迭代次数降低时间复杂度。 牛顿迭代法：其最大优点是在方程f(x) = 0的单根附近具有平方收敛，收敛速度快。而且该法还可以用来求方程的重根、复根。缺点：初值的选择会影响收敛结果。 牛顿下山法：保证函数值稳定下降，且有牛顿法的收敛速度。 有限差分方法(FDM:Finite Difference Method)是计算机数值模拟最早采用的方法，至今仍被广泛运用。该方法将求解域划分为差分网格，用有限个网格节点代替连续的求解域。有限差分法以Taylor级数展开等方法，把控制方程中的导数用网格节点上的函数值的差商代替进行离散，从而建立以网格节点上的值为未知数的代数方程组。有限差分法主要集中在依赖于时间的问题（双曲型和抛物型方程）。有限差分法方面的经典文献有Richtmeyer & Morton的《Difference Methods for Initial-Value Problems》；R. LeVeque《Finite Difference Method for Differential Equations》；《Numerical Methods for C onservation Laws》。 注：差分格式： （1）从格式的精度来划分，有一阶格式、二阶格式和高阶格式。 （2）从差分的空间形式来考虑，可分为中心格式和逆风格式。 （3）考虑时间因子的影响，差分格式还可以分为显格式、隐格式、显隐交替格式等。 目前常见的差分格式，主要是上述几种形式的组合，不同的组合构成不同的差分格式。差分方法主要适用于有结构网格，网格的步长一般根据实际地形的情况和柯朗稳定条件来决定。 构造差分的方法： 构造差分的方法有多种形式，目前主要采用的是泰勒级数展开方法。其基本的差分表达式主要有三种形式：一阶向前差分、一阶向后差分、一阶中心差分和二阶中心差分等，其中前两种格式为一阶计算精度，后两种格式为二阶计算精度。通过对时间和空间这几种不同差分格式的组合，可以组合成不同的差分计算格式。 有限差分法的不足：由于采用的是直交网格，因此较难适应区域形状的任意性，而且区分不出场函数在区域中的轻重缓急之差异，缺乏统一有效的处理自然边值条件和内边值条件的方法，难以构造高精度(指收敛阶)差分格式，除非允许差分方程联系更多的节点(这又进一步增加处理边值条件韵困难)。另外它还有编制不出通用程序的困难。 有限差分法的优点：该方法是一种直接将微分问题变为代数问题的近似数值解法，数学概念 直观，表达简单，精度可选而且在一个时间步内，对于一个给定点来说其相关的空间点只是 与该相邻的几点，而不是全部的空间点。是发展较早且比较成熟的数值方法 广义差分法（有限体积法）（GDM:Generalized Difference Method）：1953年，Mac—Neal 利用积分插值法(也称积分均衡法)建立了三角网格上的差分格 式，这就是以后通称的不规划网格上的差分法．这种方法的几何误差小，特别是给出了处理自然边值条件(及内边值条件)的有效方法，堪称差分法的一大进步。1978年，李荣华利用有限元空间和对偶单元上特征函数的推广——局部Taylor展式的公项，将积分插值法改写成广义Galerkin法形式，从而将不规则网格差分法推广为广义差分法.其基本思路是，将计算区域划分为一系列不重复的控制体积，并使每个网格点周围有 实验一：拉格朗日插值多项式 命名（源程序.cpp及工作区.dsw）：lagrange 问题： 4 //Lagrange.cpp #include 第一章 非线性方程和方程组的数值解法 1）二分法的基本原理，误差：~ 1 2 k b a x α+--< 2）迭代法收敛阶：1lim 0i p i i c εε+→∞ =≠，若1p =则要求01c << 3）单点迭代收敛定理： 定理一：若当[],x a b ∈时，[](),x a b ?∈且' ()1x l ?≤<，[],x a b ?∈，则迭代格式收敛 于唯一的根； 定理二：设()x ?满足：①[],x a b ∈时，[](),x a b ?∈， ②[]121212,,, ()(),01x x a b x x l x x l ???∈-≤-<<有 则对任意初值[]0,x a b ∈迭代收敛，且： 110 1 11i i i i i x x x l l x x x l αα+-≤ ---≤-- 定理三：设()x ?在α的邻域内具有连续的一阶导数，且'()1?α<，则迭代格式具有局部收敛性； 定理四：假设()x ?在根α的邻域内充分可导，则迭代格式1()i i x x ?+=是P 阶收敛的 () ()()0,1,,1,()0j P j P ? α?α==-≠ （Taylor 展开证明） 4）Newton 迭代法：1'() () i i i i f x x x f x +=-，平方收敛 5）Newton 迭代法收敛定理： 设()f x 在有根区间[],a b 上有二阶导数，且满足： ①：()()0f a f b <； ②：[]' ()0,,f x x a b ≠∈； ③：[]'' ,,f x a b ∈不变号 ④：初值[]0,x a b ∈使得'' ()()0f x f x <； 则Newton 迭代法收敛于根α。 第一章绪论 误差来源：模型误差、观测误差、截断误差（方法误差）、舍入误差 是的绝对误差，是的误差，为的绝对误差限（或误差限） 为的相对误差，当较小时，令 相对误差绝对值得上限称为相对误差限记为：即： 绝对误差有量纲，而相对误差无量纲 若近似值的绝对误差限为某一位上的半个单位，且该位直到的第一位非零数字共有n位，则称近似值有n位有效数字，或说精确到该位。 例：设x==…那么,则有效数字为1位，即个位上的3，或说精确到个位。 科学计数法：记有n位有效数字，精确到。 由有效数字求相对误差限：设近似值有n位有效数字，则其相对误差限为 由相对误差限求有效数字：设近似值的相对误差限为为则它有n位有效数字 令 1.x+y近似值为和的误差（限）等于误差（限）的 和 2.x-y近似值为 3.xy近似值为 4. 1．避免两相近数相减 2．避免用绝对值很小的数作除数 3．避免大数吃小数 4．尽量减少计算工作量 第二章非线性方程求根 1.逐步搜索法 设f (a) <0, f (b)> 0，有根区间为(a, b)，从x0=a出发，按某个预定步长(例如h=(b-a)/N) 一步一步向右跨，每跨一步进行一次根的搜索，即判别f(x k)=f(a+kh)的符号，若f(x k)>0(而 f(x k-1)<0),则有根区间缩小为[x k-1,x k] (若f(x k)=0，x k即为所求根), 然后从x k-1出发，把搜索步长再缩小，重复上面步骤，直到满足精度：|x k-x k-1| 列主元法 function lianzhuyuan(A,b) n=input('请输入n:') %选择阶数A=zeros(n,n); %系数矩阵A b=zeros(n,1); %矩阵b X=zeros(n,1); %解X for i=1:n for j=1:n A(i,j)=(1/(i+j-1)); %生成hilbert矩阵A end b(i,1)=sum(A(i,:)); %生成矩阵b end for i=1:n-1 j=i; top=max(abs(A(i:n,j))); %列主元 k=j; while abs(A(k,j))~=top %列主元所在行 k=k+1; end for z=1:n %交换主元所在行a1=A(i,z); A(i,z)=A(k,z); A(k,z)=a1; end a2=b(i,1); b(i,1)=b(k,1); b(k,1)=a2; for s=i+1:n %消去算法开始m=A(s,j)/A(i,j); %化简为上三角矩阵 A(s,j)=0; for p=i+1:n A(s,p)=A(s,p)-m*A(i,p); end b(s,1)=b(s,1)-m*b(i,1); end end X(n,1)=b(n,1)/A(n,n); %回代开始 for i=n-1:-1:1 s=0; %初始化s for j=i+1:n s=s+A(i,j)*X(j,1); end X(i,1)=(b(i,1)-s)/A(i,i); end X 欧拉法 clc clear % 欧拉法 p=10; %贝塔的取值 T=10; %t取值的上限 y1=1; %y1的初值 r1=1; %y2的初值 %输入步长h的值 h=input('欧拉法please input number(h=1 0.5 0.25 0.125 0.0625):h=') ; if h>1 or h<0 break end S1=0:T/h; S2=0:T/h; S3=0:T/h; S4=0:T/h; i=1; % 迭代过程 for t=0:h:T Y=(exp(-t)); R=(1/(p-1))*exp(-t)+((p-2)/(p-1))*exp(-p*t); y=y1+h*(-y1); y1=y; r=r1+h*(y1-p*r1); r1=r; S1(i)=Y; S2(i)=R; S3(i)=y; S4(i)=r; i=i+1; end t=[0:h:T]; % 红线为解析解，'x'为数值解 plot(t,S1,'r',t,S3,'x') 数值分析心得体会 篇一：学习数值分析的经验 数值分析实验的经验、感受、收获、建议班级：计算131 学号：XX014302 姓名：曾欢欢 数值分析实验主要就是学习MATLAB的使用以及对数值分析类容的应用，可以使学生更加理解和记忆数值分析学得类容，也巩固了MATLAB的学习，有利于以后这个软件我们的使用。在做实验中，我们需要具备较好的编程能力、明白MATLAB软件的使用以及掌握数值分析的思想，才能让我们独立自主的完成该作业，如果是上述能力有限的同学，需要借助MATLAB的书以及网络来完成实验。数值分析实验对于我来说还是有一定难度，所以我课下先复习了MATLAB的使用方法以及编写程序的基本类容，借助互联网和同学老师资源完成了数值分析得实验的内容。在实验书写中，我复习了各种知识，所以我认为这门课程是有必要且是有用处的，特别是需要处理大量实验数据的人员，很有必要深入了解学习它，这样在以后的工作学习里面就减少了很多计算问题也提高了实验结果的精确度。 学习数值分析的经验、感受、收获、建议数值分析的内容包括插值与逼近，数值微分与数值积分，非线性方程与线性方程组的数值解法，矩阵的特征值与特征向量计算，常微分方程数值解等。 首先我们必须明白数值分析的用途。通常所学的其他数学类学科都是由公式定理开始，从研究他们的定义，性质再到证明与应用。但实际上，尤其是工程，物理，化学等其它具体的学科。往往我们拿到 手的只是通过实验得到的数据。如果是验证性试验，需要代回到公式 进行分析，验证。但往往更多面对的是研究性或试探性试验，无具体 公式定理可代。那就必须通过插值，拟合等计算方法进行数据处理以得到一个相对可用的一般公式。还有许多计算公式理论上非常复杂，在工程中不实用，所以必须根据实际情况把它转化成多项式近似表 示。学习数值分析，不应盲目记公式，因为公事通常很长且很乏味。其次，应从公式所面临的问题以及用途出发。比如插值方法，就 是就是把实验所得的数据看成是公式的解，由这些解反推出一个近似公式，可以具有局部一般性。再比如说拟合，在插值的基础上考虑实 验误差，通过拟合能将误差尽可能缩小，之后目的也是得到一个具有 一定条件下的一般性的公式。。建议学习本门课程要结合知识与实际，比如在物理实验里面很多 数值分析的几个简单算法实现 Matlab6.5 function x=Aitken(fname,x0,N,e) k=1;e=1e-5;x2=x0+2*e; while abs(x2-x0)>e&k end disp(x); end %牛顿法 function x=newton(fname,dfname,x0,e,N) if nargin<5;N<500;end; if nargin<4;e=1e-4;end; x=x0;x0=x+2*e;k=0; while abs(x0-x)>e&k 目录 第一章非线性方程求根 (3) 1.1迭代法 (3) 1.2牛顿法 (4) 1.3弦截法 (5) 1.4二分法 (6) 第二章插值 (7) 2.1线性插值 (7) 2.2二次插值 (8) 2.3拉格朗日插值 (9) 2.4分段线性插值 (10) 2.5分段二次插值 (11) 第三章数值积分 (13) 3.1复化矩形积分法 (13) 3.2复化梯形积分法 (14) 3.3辛普森积分法 (15) 3.4变步长梯形积分法 (16) 第四章线性方程组数值法 (17) 4.1约当消去法 (17) 4.2高斯消去法 (18) 4.3三角分解法 (20) 4.4雅可比迭代法 (21) 4.5高斯—赛德尔迭代法 (23) 第五章常积分方程数值法 (25) 5.1显示欧拉公式法 (25) 5.2欧拉公式预测校正法 (26) 5.3改进欧拉公式法 (27) 5.4四阶龙格—库塔法 (28) 数值计算方法 第一章非线性方程求根 1.1迭代法 程序代码： Private Sub Command1_Click() x0 = Val(InputBox("请输入初始值x0")) ep = Val(InputBox(请输入误差限ep)) f = 0 While f = 0 X1 = (Exp(2 * x0) - x0) / 5 If Abs(X1 - x0) < ep Then Print X1 f = 1 Else x0 = X1 End If Wend End Sub 例：求f(x)=e2x-6x=0在x=0.5附近的根（ep=10-10） 1.2牛顿法 程序代码： Private Sub Command1_Click() b = Val(InputBox("请输入被开方数x0")) ep = Val(InputBox(请输入误差限ep)) f = 0 While f = 0 X1 = x0 - (x0 ^ 2 - b) / (2 * b) If Abs(X1 - x0) < ep Then Print X1 f = 1 Else x0 = X1 End If Wend End Sub 例：求56的值。（ep=10-10） 上机实习题一 一、题目： 已知A 与b 12.38412,2.115237,-1.061074,1.112336,-0.113584,0.718719,1.742382,3.067813,-2.031743 2.115237,19.141823,- 3.125432,-1.012345,2.189736,1.563849,-0.784165,1.112348,3.123124 -1.061074,-3.125432,15.567914,3.123848,2.031454,1.836742,-1.056781,0.336993,-1.010103 1.112336,-1.012345,3.123848,27.108437, 4.101011,-3.741856,2.101023,-0.71828,-0.037585 A= -0.113584,2.189736,2.031454,4.101011,19.897918,0.431637,-3.111223,2.121314,1.784137 0.718719,1.563849,1.836742,-3.741856,0.431637,9.789365,-0.103458,-1.103456,0.238417 1.742382,-0.784165,-1.056781,2.101023,-3.111223,-0.103458,14.7138465,3.123789,-2.213474 3.067813,1.112348,0.336993,-0.71828,2.121314,-1.103456,3.123789,30.719334,4.446782 -2.031743,3.123124,-1.010103,-0.037585,1.784317,0.238417,-2.213474,4.446782,40.00001 b={2.1874369,33.992318,-25.173417,0.84671695,1.784317,-86.612343,1.1101230,4.719345,-5.6784392} 1．用Household 变换，把A 化为三对角阵B （并打印B ）。 2．用超松弛法求解BX=b （取松弛因子ω=1.4,X (0)=0,迭代9次）。 3．用列主元素消去法求解BX=b 。 二、解题方法的理论依据: 1 、用Householder 变换的理论依据 ﹝1﹞ 令A0=A,a(ij)1=a(ij),已知Ar_1即Ar_1=a(ij)r ﹝2﹞ Sr=sqrt(pow(a,2)) ﹝3﹞ a(r)=Sr*Sr+abs(a(r+1,r))*Sr ﹝4﹞ y(r)=A(r_1)*u?/a? ﹝5﹞ Kr=(/2)*Ur 的转置*Yr/a? ﹝6﹞ Qr=Yr-Kr*Ur ﹝7﹞ Ar=A(r-1)-(Qr*Ur 的转置+Ur*Qr 的转置) r=1,2,,……,n -2 2 、用超松弛法求解 其基本思想：在高斯方法已求出x (m),x (m-1) 的基础上，组合新的序列，从而加快收敛速度。 其算式： ??? ??? ?=+?=+-+----=][][0][0][~ ][]][[][]1[0]][[]1][[]1[]][[]1][[][~ i X i X i X i X w i X i i B i b i X i i B i i B i X i i B i i B i X 其中ω是超松弛因子，当ω>1时，可以加快收敛速度 3 、用消去法求解 用追赶消去法求Bx=b 的方法： ][]1[1i b i d =+ , ]][1[]2[1i i B i a +=+ , ][]1[]1][[][]][[]1[]1][[i b i X i i B i X i i B i X i i B =+?++?+-?- 《计算方法》实验内容 一．实验一：用两种不同的顺序计算 644834.110000 1 2 ≈∑=-n n ，分析其误差的变化。 1.实验目的：通过正序反序两种不同的顺序求和，比较不同算法的误差；了解在计算机中大数吃小数的现象，以后尽量避免；体会单精度和双精度数据的差别。 2.算法描述：累加和s=0； 正序求和： 对于n=1,2,3,......,10000 s+=1.0/(n*n); 反序求和： 对于n=10000,9999,9998,.....,1 s+=1.0/(n*n); 3.源程序： #双精度型# #include 4.运行结果： 双精度型运行结果： 单精度型运行结果： 5.对算法的理解与分析：舍入误差在计算机中会引起熟知的不稳定，算法不同，肯结果也会不同，因此选取稳定的算法很重要。选取双精度型数据正反序求和时结果一致，但选用单精度型数据时，求和结果不一致，明显正序求和结果有误差，所以第一个算法较为稳定可靠。 二．实验二： 1、拉格朗日插值 按下列数据 x i -3.0 -1.0 1.0 2.0 3.0 y i 1.0 1.5 2.0 2.0 1.0 作二次插值，并求x 1=-2，x 2 =0，x 3 =2.75时的函数近似值 2牛顿插值 按下列数据 x i 0.30 0.42 0.50 0.58 0.66 0.72 y i 1.0440 3 1.0846 2 1.1180 3 1.1560 3 1.19817 1.23223 作五次插值，并求x 1=0.46，x 2 =0.55，x 3 =0.60时的函数近似值. 1.实验目的：通过拉格朗日插值和牛顿插值的实例，了解两种求解方法，并分析各自的优缺点。 2.算法描述： 3.源程序： 拉格朗日插值： #include 《数值分析》实验报告 第二章=解线性方程组的直接方法 2、试用MATLAB软件编程实现追赶法求解三对角方程组的算法, 并考虑梯形电阻电路问题，电路如下： 其中电路中的各个电流f2?…，珀｝须满足下列线性方程组: 一耳十务一巧=0 一巧4■笃一為=0 一迢+呂一逼=0 一耳+兀_%=0 一爲+5% 一巧=0 -2f6+5f7-2f8 =0 -2i7+5r g=0 设F = 220V,氏= 27Q,运用求各段电路的电流量。 解：v!R = — ?8_1481 27 上述方程组可用矩阵表示为： ?2_2000000h81481 _25_200000h0 0_25-20000h0 00_25_2000h0 000-25_200h0 0000_25-20S0 00000-25_2h0 000000-25?8.0 認翟蚁蛊 X -p u a) A .H )q /(i p n .H )p T (.rH )p T 二—二丄岂4 晒口?〈狠* A8)q 、 (8)PH(8)p -p u a) -U —.H ) P *-H ) F (.H ) P A .H ) P 丄) y (.H )q A .H )q Q 丄)q/u)y(2 晒口§前厂< 00 勺：。 4 坦祖咽鱼磋￡紙夜 lll < =0 0 0 0 0 0 0 T 畧 oo v p 二甲 d — d — 甲甲甲 d — v 。 二g g g g g g g d v q 氏 —甲 甲甲甲甲 d — O V B 曲恥 qedelN 寸 61—<10905 30p ^ s ￡0- S S I 黑理q 。脅 u) 旺=味蛋二 第一章绪论（1-4） 一、误差来源及分类 二、误差的基本概念 1.绝对误差及绝对误差限 2.相对误差及相对误差限 3.有效数字 三、数值计算的误差估计 1.函数值的误差估计 2.四则运算的误差估计 四、数值计算的误差分析原则 第二章插值（1.2.4-8） 一、插值问题的提法（定义）、插值条件、插值多项式的存在唯一性 二、拉格朗日插值 1.拉格朗日插值基函数的定义、性质 2.用拉格朗日基函数求拉格朗日多项式 3.拉格朗日插值余项（误差估计） 三、牛顿插值 1.插商的定义、性质 2.插商表的计算 3.学会用插商求牛顿插值多项式 四、等距节点的牛顿插值 1.差分定义、性质及计算（向前、向后和中心） 2.学会用差分求等距节点下的牛顿插值公式 五、学会求低次的hermite插值多项式 六、分段插值 1.分段线性插值 2.分段三次hermite插值 3.样条插值 第三章函数逼近与计算（1-6） 一、函数逼近与计算的提法（定义）、常用两种度量标准（一范数、二范数\平方逼近） 二、基本概念 连续函数空间、最佳一次逼近、最佳平方逼近、内积、内积空间、偏差与最小偏差、偏差点、交错点值、平方误差 三、学会用chebyshev定理求一次最佳一致逼近多项式，并估计误差（最大偏差） 四、学会在给定子空间上通过解方程组求最佳平方逼近，并估计误差（平方误差） 五、正交多项式（两种）定义、性质，并学会用chebyshev多项式性质求特殊函数的（降阶）最佳一次逼近多项式 六、函数按正交多项式展开求最佳平方逼近多项式，并估计误差 七、一般最小二乘法（多项式拟合）求线性拟合问题 第四章数值分析（1-4） 一、数值求积的基本思想及其机械求积公式 《计算方法》期中复习试题 一、填空题： 1、已知3.1)3(,2.1)2(,0.1)1(===f f f ，则用辛普生（辛卜生）公式计算求得 ?≈3 1 _________ )(dx x f ,用三点式求得≈')1(f 。 答案：2.367，0.25 2、1)3(,2)2(,1)1(==-=f f f ，则过这三点的二次插值多项式中2 x 的系数为 ，拉 格朗日插值多项式为 。 答案：-1， )2)(1(21 )3)(1(2)3)(2(21)(2--------= x x x x x x x L 3、近似值*0.231x =关于真值229.0=x 有( 2 )位有效数字； 4、设)(x f 可微,求方程)(x f x =的牛顿迭代格式是( )； 答案 )(1)(1n n n n n x f x f x x x '--- =+ 5、对1)(3 ++=x x x f ,差商=]3,2,1,0[f ( 1 ),=]4,3,2,1,0[f ( 0 )； 6、计算方法主要研究( 截断 )误差和( 舍入 )误差； 7、用二分法求非线性方程 f (x )=0在区间(a ,b )内的根时，二分n 次后的误差限为 ( 1 2+-n a b )； 8、已知f (1)＝2，f (2)＝3，f (4)＝5.9，则二次Newton 插值多项式中x 2系数为( 0.15 )； 11、 两点式高斯型求积公式?1 d )(x x f ≈( ?++-≈1 )] 321 3()3213([21d )(f f x x f )，代数精度 为( 5 )； 12、 为了使计算 32)1(6 )1(41310-- -+-+ =x x x y 的乘除法次数尽量地少，应将该表达 式改写为 11 ,))64(3(10-= -++=x t t t t y ，为了减少舍入误差，应将表达式1999 2001- 第一章 引论（习题） 2．证明：x 的相对误差约等于x 的相对误差的1/2. 证明 记 x x f = )( ，则 ) ()(* ** x x x x x x x x f E r +-= -= )(21**x E x x x x x x r ≈-?+= . □ 3．设实数a 的t 位β进制浮点机器数表示为)(a fl . 试证明 t b a b a fl -≤ +*=*12 1||),1/()()(βδδ， 其中的记号*表示+、-、?、/ 中一种运算. 证明： 令： ) () ()(b a fl b a fl b a **-*= δ 可估计： 1|)(|-≥*c b a fl β （c 为b a *阶码）， 故： 121||--≤ c t c ββδt -=12 1β 于是： )1()()(δ+*=*b a b a fl . □ 4．改变下列表达式使计算结果比较精确： （1） ;1||, 11211<<+--+x x x x 对 （2） ;1,11>>- -+ x x x x x 对 （3） 1||,0,cos 1<<≠-x x x x 对. 解 (1) )21()1(22 x x x ++. (2) ) 11(2x x x x x -++. (3) x x x x x x x cos 1sin )cos 1(sin cos 12+≈+=-. □ 6．设937.0=a 关于精确数x 有3位有效数字，估计a 的相对误差. 对于x x f -=1)(，估计)(a f 对于)(x f 的误差和相对误差. 解 a 的相对误差：由于 31021|)(|-?≤ -=a x x E . x a x x E r -=)(, 221018 1 10921)(--?=?≤ x E r . （1Th ） )(a f 对于)(x f 的误差和相对误差. |11||)(|a x f E ---== ()25 .0210 11321??≤ -+---a x x a =3 10- 33 104110 |)(|--?=-≤a f E r . □ 9．序列}{n y 满足递推关系：1101.100-+-=n n n y y y . 取01.0,110 ==y y 及 01.0, 101150=+=-y y ，试分别计算5y ，从而说明该递推公式对于计算是不稳 定的. 解 递推关系： 1101.100-+-=n n n y y y (1) 取初值 10=y ， 01.01=y 计算 可得： 110 01.1002 2-?=-y 10001.1-=410-= 6 310-=y ， 8 410 -=y ， 10 510-=y ， … (2) 取初值 5 0101-+=y ， 2 110 -=y , 记： n n n y y -=ε, 序列 {}n ε ，满足递推关系，且 5 010--=ε ， 01=ε 1101.100-+-=n n n εεε, 于是： 5210-=ε, 531001.100-?=ε, 55241010)01.100(---?=ε, 5 5351002.20010)01.100(--?-?=ε, 变步长梯形求积法计算定积分 1．原理： 变步长求积法的主要思想是利用若干小梯形的面积代替原方程的积分，当精度达不到要求时，可以通过增加点数对已有的区间再次划分，达到所需精度时即可；其中由于新的式子中有原来n点中的部分项，故可以省略一些计算，符合了计算机计算存储的思想。 主要公式：T2n=T n/2+(h/2)*Σf(x k+; 2．C++语言实现方式： 通过每次的T n值和新增的函数值点计算T2n，再通过判断|T n-T2n|的大小来判断是否达到精度要求。 3．源程序如下： #include"" #include"" double f(double x)//预先输入的待积分函数 { double s; s=log(x*x); return(s); } double ffts(double a,double b,double eps) { int n,k; double fa,fb,h,t1,p,s,x,t; fa=f(a); fb=f(b); n=1; h=b-a; t1=h*(fa+fb)/2; p=eps+1; while(p>=eps) { s=0; for(k=0;k<=n-1;k++) { x=a+(k+*h; s=s+f(x); } t=t1/2+h*s/2; p=fabs(t1-t); cout<<"步长n为:"< 计算流体力学常用数值方法简介 李志印 熊小辉 吴家鸣 (华南理工大学交通学院) 关键词 计算流体力学 数值计算 一 前 言 任何流体运动的动力学特征都是由质量守恒、动量守恒和能量守恒定律所确定的,这些基本定律可以由流体流动的控制方程组来描述。利用数值方法通过计算机求解描述流体运动的控制方程,揭示流体运动的物理规律,研究流体运动的时一空物理特征,这样的学科称为计算流体力学。 计算流体力学是一门由多领域交叉而形成的一门应用基础学科,它涉及流体力学理论、计算机技术、偏微分方程的数学理论、数值方法等学科。一般认为计算流体力学是从20世纪60年代中后期逐步发展起来的,大致经历了四个发展阶段:无粘性线性、无粘性非线性、雷诺平均的N-S方程以及完全的N-S方程。随着计算机技术、网络技术、计算方法和后处理技术的迅速发展,利用计算流体力学解决流动问题的能力越来越高,现在许多复杂的流动问题可以通过数值计算手段进行分析并给出相应的结果。 经过40年来的发展,计算流体力学己经成为一种有力的数值实验与设计手段,在许多工业领域如航天航空、汽车、船舶等部门解决了大量的工程设计实际问题,其中在航天航空领域所取得的成绩尤为显著。现在人们已经可以利用计算流体力学方法来设计飞机的外形,确定其气动载荷,从而有效地提高了设计效率,减少了风洞试验次数,大大地降低了设计成本。此外,计算流体力学也己经大量应用于大气、生态环境、车辆工程、船舶工程、传热以及工业中的化学反应等各个领域,显示了计算流体力学强大的生命力。 随着计算机技术的发展和所需要解决的工程问题的复杂性的增加,计算流体力学也己经发展成为以数值手段求解流体力学物理模型、分析其流动机理为主线,包括计算机技术、计算方法、网格技术和可视化后处理技术等多种技术的综合体。目前计算流体力学主要向二个方向发展:一方面是研究流动非定常稳定性以及湍流流动机理,开展高精度、高分辩率的计算方法和并行算法等的流动机理与算法研究;另一方面是将计算流体力学直接应用于模拟各种实际流动,解决工业生产中的各种问题。 二 计算流体力学常用数值方法 流体力学数值方法有很多种,其数学原理各不相同,但有二点是所有方法都具备的,即离散化和代数化。总的来说其基本思想是:将原来连续的求解区域划分成网格或单元子区 (1) 解题过程如下： （1）MATLAB中创建复化梯形公式和复化辛普森公式的 M 文件：1）复化梯形公式文件： function s=T_fuhua(f,a,b,n) h=(b-a)/n; s=0; for k=1:(n-1) x=a+h*k; s=s+feval(f,x); end s=h*(feval(f,a)+feval(f,b))/2+h*s; 2）复化辛普森公式文件： function s=S_fuhua(f,a,b,n) h=0; h=(b-a)./(2*n); s1=0; https://www.doczj.com/doc/055110101.html, -5- s2=0; for k=1:n-1 x=a+h*2*k; s1=s1+feval(f,x); end for k=1:n x=a+h*(2*k-1); s2=s2+feval(f,x); end s=h*(feval(f,a)+feval(f,b)+s1*2+s2*4)/3; 在MATLAB中输入： f=inline('x/(4+x^2)');a=0;b=1; %inline 构造内联函数对象 for n=2:10 s(n-1)=T_fuhua(f,a,b,n);s(n-1)=vpa(s(n-1),10); %调用复化梯形公式，生成任意精度的数值 end exact=int('x/(4+x^2)',0,1);exact=vpa(exact,10) %求出积分的精确值 输出结果：exact = .1115717755 s = Columns 1 through 6 0.1088 0.1104 0.1109 0.1111 0.1113 0.1114 Columns 7 through 9 0.1114 0.1114 0.1115 在MATLAB中输入以下函数用以画出计算误差与 n 之间的曲线： r=abs(exact-s); n=2:10; plot(double(n),double(r(n-1))) 得到结果如图所示： （2）在 MATLAB中输入以下程序代码： f=inline('x/(4+x^2)');a=0;b=1;n=9; %inline 构造内联函数对象 t=T_fuhua(f,a,b,n);t=vpa(t,10) s=S_fuhua(f,a,b,n);s=vpa(s,10)数值分析综述-《数值分析与算法》徐士良
数值计算方法比较
数值分析所有代码
数值分析重点公式
数值分析(计算方法)总结
数值分析的MATLAB程序
数值分析心得体会
【免费下载】数值分析的几个简单算法实现
数值计算方法大作业
数值分析编程题c语言
数值分析计算方法
数值分析——编程作业..docx
数值分析知识点
数值计算方法》试题集及答案
数值分析习题汇总
数值分析与算法变步长梯形求积法计算定积分
计算流体力学常用数值方法简介[1]
数值分析 matlab 实验4
相关主题
文本预览