所有的算法 数值方法.docx
- 文档编号:18022767
- 上传时间:2023-08-05
- 格式:DOCX
- 页数:33
- 大小:115.17KB
所有的算法 数值方法.docx
《所有的算法 数值方法.docx》由会员分享,可在线阅读,更多相关《所有的算法 数值方法.docx(33页珍藏版)》请在冰点文库上搜索。
所有的算法数值方法
数值分析上机实验
一、三次样条插值:
P52
第一种边界条件:
9、给定函数
的函数表和边界条件
,求三次样条插值函数
,并求
的近似值。
函数表
75
76
77
78
79
80
2.768
2.833
2.903
2.979
3.062
3.153
源代码:
#include
#include
intmain()
{constintn=6;//给定节点个数
doublex[n]={75,76,77,78,79,80};
doubley[n]={2.768,2.833,2.903,2.979,3.062,3.153};//初始化插值节点xi、yi
doublexx=78.3;//点xx
doubleAlph[n],Bet[n];//α、β
doublea[n],b[n];
doublesx;//三次样条插值函数s(x)在点xx处的函数值
//计算h[i]=x[i+1]+x[i](i=0,1,...,n-2)
doubleh[n-1];
for(inti=0;i { h[i]=x[i+1]-x[i]; } //计算αi和βi(i=0,1,...,n-1),对应第二种边界条件 Alph[0]=1; Alph[n-1]=0; Bet[0]=3/h[0]*(y[1]-y[0]); Bet[n-1]=3/h[n-2]*(y[n-1]-y[n-2]); for(i=1;i { Alph[i]=h[i-1]/(h[i-1]+h[i]); Bet[i]=3*((1-Alph[i])/h[i-1]*(y[i]-y[i-1])+Alph[i]/h[i]*(y[i+1]-y[i])); } //计算ai和bi(i=0,1,...,n-1) a[0]=-Alph[0]/2;b[0]=Bet[0]/2; for(i=1;i { a[i]=-Alph[i]/(2+(1-Alph[i])*a[i-1]); b[i]=(Bet[i]-(1-Alph[i])*b[i-1])/(2+(1-Alph[i])*a[i-1]); } //计算mi; doublem[n+1]; m[n]=0; for(i=n-1;i>=0;i--) { m[i]=a[i]*m[i+1]+b[i]; } //判断xx所在的区间[xi,xi+1], i=0; while (1) { if(xx>=x[i]&&xx break; elsei++; } //求出s(x)在点xx处的值sx并输出 sx=(1+2*(xx-x[i])/(x[i+1]-x[i]))*pow((xx-x[i+1])/(x[i]-x[i+1]),2)*y[i] +(1+2*(xx-x[i+1])/(x[i]-x[i+1]))*pow((xx-x[i])/(x[i+1]-x[i]),2)*y[i+1] +(xx-x[i])*pow((xx-x[i+1])/(x[i]-x[i+1]),2)*m[i] +(xx-x[i+1])*pow((xx-x[i])/(x[i+1]-x[i]),2)*m[i+1]; cout<<"函数表\n" <<"x757677787980\n" <<"y=f(x)2.7682.8332.9032.9793.0623.153\n" <<"s\"(75)=0,s\"(80)=0\n" <<"由三次样条插值的第二种边界条件可得f(78.3)的近似值为: "< return0; } 运行结果截图: 第二种边界条件: 10、给定函数 的函数表和边界条件 ,求三次样条插值函数 ,并求 的近似值。 函数表 0.25 0.3 0.39 0.45 0.53 0.5 0.5477 0.6245 0.6708 0.728 源代码: #include #include intmain() {constintn=5;//给定节点个数 doublex[n]={0.25,0.3,0.39,0.45,0.53}; doubley[n]={0.5,0.5477,0.6245,0.6708,0.728};//初始化插值节点xi、yi doublexx=0.35;//点xx doublem0=1,mn=0.6868;//s'(x0)=m0,s'(xn)=mn doubleAlph[n],Bet[n];//α、β doublea[n],b[n]; doublesx;//三次样条插值函数s(x)在点xx处得函数值 //计算h[i]=x[i+1]+x[i](i=0,1,...,n-2) doubleh[n-1]; for(inti=0;i { h[i]=x[i+1]-x[i]; } //计算αi和βi(i=0,1,...,n-1),对应第一种边界条件 Alph[0]=0; Alph[n-1]=1; Bet[0]=2*m0; Bet[n-1]=2*mn; for(i=1;i { Alph[i]=h[i-1]/(h[i-1]+h[i]); Bet[i]=3*((1-Alph[i])/h[i-1]*(y[i]-y[i-1]) +Alph[i]/h[i]*(y[i+1]-y[i])); } //计算ai和bi(i=0,1,...,n-1) a[0]=-Alph[0]/2;b[0]=Bet[0]/2; for(i=1;i { a[i]=-Alph[i]/(2+(1-Alph[i])*a[i-1]); b[i]=(Bet[i]-(1-Alph[i])*b[i-1])/(2+(1-Alph[i])*a[i-1]); } //计算mi(i=0,1,...,n) doublem[n+1]; m[n]=0; for(;i>=0;i--) { m[i]=a[i]*m[i+1]+b[i]; } //判断xx所在的区间[xi,xi+1] i=0; while (1) { if(xx>=x[i]&&xx break; elsei++; } //求出s(x)在点xx处的值sx并输出 sx=(1+2*(xx-x[i])/(x[i+1]-x[i]))*pow((xx-x[i+1])/(x[i]-x[i+1]),2)*y[i] +(1+2*(xx-x[i+1])/(x[i]-x[i+1]))*pow((xx-x[i])/(x[i+1]-x[i]),2)*y[i+1] +(xx-x[i])*pow((xx-x[i+1])/(x[i]-x[i+1]),2)*m[i] +(xx-x[i+1])*pow((xx-x[i])/(x[i+1]-x[i]),2)*m[i+1]; cout<<"函数表\n" <<"x0.250.30.390.450.53\n" <<"y=f(x)0.50.54770.62450.67080.728\n" <<"s\'(0.25)=1,s\'(0.53)=0.6868\n" <<"由三次样条插值的第一种边界条件可得f(0.35)的近似值为: "< return0; } 运行结果截图: 二、自动选取步长梯形法: P97 7、使用自动选取步长梯形法计算积分 的近似值。 (给定ε=0.01) 源代码: #include #include doublefun_t(doublet)//被积函数 { return2/(1+t*t); } intmain() { doublea,b,Eps; cout<<"请输入积分上下限a,b及精度Eps: "; cin>>a>>b>>Eps; doubleh=(b-a)/2; intn=1; doubleT0,T1=(fun_t(a)+fun_t(b))*h;//T0为前次积分近似值,T1为后次积分近似值 doubleS; do { T0=T1; S=0; for(intk=1;k<=n;k++) { S=S+fun_t(a+(2*k-1)*h/n); } T1=T0/2+S*h/n; n=2*n; }while(fabs(T1-T0)>=3*Eps); cout<<"所求积分的近似值为: "< return0; } 运行结果截图: 三、Romberg求积法: P97 8、使用Romberg求积法计算积分 的近似值。 (给定ε=0.01,且取 ) 源代码: #include #include staticdoubleT[200][200]; doublefun_x(doublex)//被积函数 { returnsqrt(x); } doubleRomberg(doublea,doubleb,doubleEps) { intk=0;//用来记录把区间[a,b]2等分的次数 T[0][0]=(b-a)/2*(fun_x(a)+fun_x(b)); do { k++; doubletemp=0; for(inti=1;i<=pow(2,k-1);i++) { temp+=fun_x(a+(2*i-1)*(b-a)/pow(2,k)); } T[0][k]=0.5*(T[0][k-1]+((b-a)/pow(2,k-1))*temp); for(intm=1;m<=k;m++) { T[m][k-m]=(pow(4,m)*T[m-1][k-m+1]-T[m-1][k-m])/(pow(4,m)-1); } }while(fabs(T[k][0]-T[k-1][0])>=Eps); returnT[k][0]; } intmain() { doublea,b,Eps; cout<<"请输入积分上下限a,b及精度Eps: "; cin>>a>>b>>Eps; cout<<"所求积分的近似值为: "< return0; } 运行结果截图: 四、列主元高斯消去法: 所求线性方程组为: = 溢出控制: ε=0.0001 源代码: #include #include conststaticintN=4; staticdoubleA[N][N]={2,-1,0,0,-1,2,-1,0,0,-1,2,-1,0,0,-1,2};//所解线性方程组的系数矩阵 staticdoubleb[N]={1,0,1,0};//所解线性方程组的右端项 staticdoublex[N];//未知数向量 doubleEps=0.0001;//溢出控制 voidswap(double&i,double&j)//交换数据 { doubletemp=i; i=j; j=temp; } voidprint(double*x,intn)//打印 { for(inti=0;i cout<<"x"< } voidmain() { intmax=0;//存放主元行下标 doubleT=0;//存放主元 for(intk=0;k {//选主元 T=fabs(A[k][k]); max=k; for(inti=k+1;i { if(T { T=fabs(A[i][k]); max=i; } } if(T { cout<<"求解失败! "< return; } else { if(max! =k) { doubletemp; for(inti=0;i { temp=A[max][i]; A[max][i]=A[k][i]; A[k][i]=temp; } swap(b[max],b[k]); } } //消元 for(i=k+1;i { T=A[i][k]/A[k][k]; b[i]=b[i]-T*b[k]; for(intj=k;j A[i][j]=A[i][j]-T*A[k][j]; } } //回代 if(fabs(A[N-1][N-1])<=Eps) { cout<<"求解失败! "< return; } else { x[N-1]=b[N-1]/A[N-1][N-1]; for(inti=N-2;i>=0;i--) { doubleT=0; for(intj=i+1;j T+=A[i][j]*x[j]; x[i]=(b[i]-T)/A[i][i]; } } print(x,N); } 实验结果截图: 五、列主元LU消去法 所求线性方程组为: = 溢出控制: ε=0.0001 源代码: include #include conststaticintN=4; staticdoubleA[N][N]={2,-1,0,0,-1,2,-1,0,0,-1,2,-1,0,0,-1,2};//线性方程组系数矩阵 staticdoubleb[N]={1,0,1,0};//右端项向量 staticdoubleY[N]; staticdoubleX[N];//未知数向量 staticdoubleS[N];//用于选取列主元 doubleEps=0.0001; voidmain() { for(intk=0;k {//选列主元 intindex=k; for(inti=k;i { doubletemp=0; for(intm=0;m { temp=temp+A[i][m]*A[m][k]; } S[i]=A[i][k]-temp; if(S[index] { index=i; } } if(fabs(S[index]) { cout<<"求解失败! "< return; } else {//交换行 doubletemp; for(i=k;i { temp=A[index][i]; A[index][i]=A[k][i]; A[k][i]=temp; } temp=b[index]; b[index]=b[k]; b[k]=temp; } //构造L、U矩阵 for(intj=k;j { doubletemp=0; for(intm=0;m { temp=temp+A[k][m]*A[m][j]; } A[k][j]=A[k][j]-temp;//先构造U一行的向量 } for(i=k+1;i { doubletemp=0; for(intm=0;m { temp=temp+A[i][m]*A[m][k]; } A[i][k]=(A[i][k]-temp)/A[k][k];//再构造L一列的向量 } } //求解LY=b Y[0]=b[0]; for(inti=1;i { doubletemp=0; for(intj=0;j { temp=temp+A[i][j]*Y[j]; } Y[i]=b[i]-temp; } //求解UX=Y X[N-1]=Y[N-1]/A[N-1][N-1]; for(i=N-2;i>=0;i--) { doubletemp=0; for(intj=i+1;j { temp=temp+A[i][j]*X[j]; } X[i]=(Y[i]-temp)/A[i][i]; } //打印X cout<<"该线性方程组的解为: "< for(i=0;i { cout<<"x"< } } 源代码: 六、简单迭代法(Jacobi迭代) 所求线性方程组为: 初始向量为: = Y= 容许误差: ε=0.0001 源代码: #include #include conststaticintN=4;//矩阵维数 staticdoubleA[N][N]={2,-1,0,0,-1,2,-1,0,0,-1,2,-1,0,0,-1,2};//所解线性方程组的系数矩阵; staticdoubleb[N]={1,0,1,0};//所解线性方程组的右端项; staticdoubleX[N]={0};//未知数向量 staticdoubleY[N]={0};//初始向量 staticdoubleg[N]={0}; doubleEps=0.0001;//容许误差 intNorm(double*X,double*Y)//判断矩阵X-Y的范数的精度,满足精度返回1,否则返回0 { doubletemp=0; for(inti=0;i temp=temp+fabs(X[i]-Y[i]); if(temp return1; else return0; } voidmain() { intk=1;//迭代次数 intM=50;//允许最大迭代次数 doubleT; for(inti=0;i {//形成迭代矩阵存放在A中 if(fabs(A[i][i]) { cout<<"求解失败"< return; } else T=A[i][i]; for(intj=0;j { A[i][j]=-A[i][j]/T; A[i][i]=0; g[i]=b[i]/T; } } do{//迭代 for(inti=0;i { doubletemp=0; for(intj=0;j { temp=temp+A[i][j]*Y[j]; } X[i]=g[i]+temp; } if(Norm(X,Y)==1) { cout<<"该线性方程组的解为: "< for(i=0;i { cout<<"X"< } cout<<"迭代次数为: "< cout<<"容许误差为: "< cout<<"允许最大迭代次数为: "< break; } else { if(k { k++; for(inti=0;i Y[i]=X[i]; } else { cout<<"求解失败! "< break; } } }while (1); return; } 实验结果截图: 1.当迭代次数在允许迭代次数以内: 2.当迭代次数超出允许迭代次数: 七、Seidel迭代 所求线性方程组为: 初始向量为: = Y= 容许误差: ε=0.0001 源代码: #include #include conststaticintN=4;//矩阵维数 stat
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 所有的算法 数值方法 所有 算法 数值 方法
![提示](https://static.bingdoc.com/images/bang_tan.gif)