数值计算课程设计(陕西科技大学)

1、经典四阶龙格库塔法解一阶微分方程1.1、 算法说明龙格-库塔（Runge-Kutta）方法是一种在工程上应用广泛的高精度单步算法。由于此算法精度高，米取措施对误差进行抑制，所以其实现原理也较复杂。该算法是构建在数学支持的基础之上的。4阶龙格-库塔方法（RK4）可模拟N=4的泰勒方法的精度。这种算法可以描述为,自初始点5yo）开始，利用下面的计算方法生成近似序列『 h儿+im石血十曲十2怠十短）%=扛耳’闰F%十卜必十+加1） （］_］）怠=/（输十扣兀十期』%=珥十徐兀十曲心）1.2、 经典四阶龙格库塔法解一阶微分方程算法流程图图1-1经典四阶龙格库塔法解一阶微分方程算法流程图

1.3、经典四阶龙格库塔法解一阶微分方程程序调试Ii输人初值斓，购沖2牖人区囘:13i输入步彳矣：0.1|626.11.904917708Q.21.819349692Q.31.742831934Q.41.674923048Q.51.615203175Q.61.5632729310.71.5187524120.S1.4812802420.91.45051266911.4261227051.11.407799311.21.395246616i1.31.3881831811.41.3863412911.5i.3894£628S1.£1.3973159351.71.4096598971.61.4262787181.91.44696417421・471517845厂i人取_ytocontinue图1-2经典四阶龙格库塔法解一阶微分方程程序调试1.4、经典四阶龙格库塔法解一阶微分方程代码#include<iostream>#include<iomanip>usingnamespacestd;//f为函数的入口地址，xO、yO为初值，xn为所求点，step为计算次数doubleRunge_Kuta(double(*f)(doublex,doubley),doublex0,doubley0,doublexn,intstep){doublekl,k2,k3,k4,result;doubleh=(xn-x0)/step;if(step<=0)return(y0);if(step==1){k1=f(x0,y0);k2=f(x0+h/2,y0+h*k1/2);k3=f(x0+h/2,y0+h*k2/2);k4=f(x0+h,y0+h*k3);result=y0+h*(k1+2*k2+2*k3+k4)/6;}else{doublex1,y1;x1=xn-h;y1=Runge_Kuta(f,x0,y0,xn-h,step-1);k1=f(x1,y1);k2=f(x1+h/2,y1+h*k1/2);k3=f(x1+h/2,y1+h*k2/2);k4=f(x1+h,y1+h*k3);result=y1+h*(k1+2*k2+2*k3+k4)/6;}return(result);}intmain(){doublef(doublex,doubley);doublex0,y0;doublea,b;//intstep;cout<<"请输入初值xO,yO:";cin>>x0>>y0;cout<<"请输入区间：";cin>>a>>b;//doublex0=0,y0=1;doublex,y,step;inti;cout<〈“请输入步长：";cin>>step;//step=0.1;cout.precision(10);for(i=0;i<=(b-a)/step;i++){x=x0+i*step;cout<<setw(8)<<x<<setw(18)<<Runge_Kuta(f,x0,y0,x,i)<<endl;}return(0);}doublef(doublex,doubley){doubler;r=(x-y)/2;return(r);1010std;A[m][m],b[m],Aug[m][m+1];b[m],doublex[m],intn)i,i1,j,k;Aug[m][m+1],maxele,Temp,l,s;//构建增广矩阵Aug(j=0;j<n;j++)//输出增广矩阵Aug(j=0;j<n;j++)〃.7//求主元2、高斯列主元法解线性方程组2.1、 算法说明首先将线性方程组做成增光矩阵，对增广矩阵进行行变换。对第元素aii，在第i列中，第i行及以下的元素选取绝对值最大的元素，将该元素最大的行与第i行交换，然后采用高斯消元法将新得到的aii消去第i行以下的元素。一次进行直到 ann。从而得到上三角矩阵。再对得到的上三角矩阵进行回代操作，即可以得到方程组的解。2.2、 高斯列主元算法代码#include<iostream>#include<iomanip>#include<cmath>#define musing namespacedoublevoidequ(doubleA[m][m],double{intdoublefor(i=0;i<n;i++){forAug[i][j]=A[i][j];Aug[i][n]=b[i];}for(i=0;i<n;i++){forcout<<Aug[i][j]<<"cout<<Aug[i][n]<<endl;}for(i1=0;i1<n-1;i1++){maxele=fabs(Aug[i1][i1]);k=i1;for(i=i1;i<n;i++)if{maxele=fabs(Aug[i][i1]);k=i;cout<<"i1="<<i1<<""<<"maxele="<<maxele<<"}if{maxele=fabs(Aug[i][i1]);k=i;cout<<"i1="<<i1<<""<<"maxele="<<maxele<<"}(maxele<fabs(Aug[i][i1]))for(j=i1;j<n+1;j++){Temp=Aug[i1][j];Aug[i1][j]=Aug[k][j];Aug[k][j]=Temp;}for(i=0;i<n;i++){forcout<<Aug[i][j]<<"cout<<Aug[i][n]<<endl;}"<<"k="<<k<<""<<Aug[i][i1]<<""<<endl;//换行for{l=-Aug[k][i1]/Aug[i1][i1];forAug[k][j]=Aug[k][j]+l*Aug[i1][j];}cout<<"每次消元结束for(i=0;i<n;i++){forcout<<Aug[i][j]<<"cout<<Aug[i][n]<<endl;//换行后输出增广矩阵Aug(j=0;j<n;j++)〃.7(k=i1+1;k<n;k++)(j=i1;j<n+1;j++)输出增广矩阵Aug"<<endl;//每次消元结束输出增广矩阵Aug(j=0;j<n;j++)x[n-1]=Aug[n-1][n]/Aug[n-1][n-1];(i=n-2;i>=0;i=i-1)(j=i+1;j<n;j++)int{voidequ(doubledoubleintcout<<"(i=n-2;i>=0;i=i-1)(j=i+1;j<n;j++)int{voidequ(doubledoubleintcout<<"输cin>>n;if{cout<<"return}for{cout<<"请输for(j=0;j<n;j++)cin>>A[i][j];cout<<"请cin>>b[i];}入问题规模太大A[m][m],double未知量的需更改原程序0入A的第输入b的第b[m],doublemain()x[m],intn);A[m][m],b[m],x[m];i,j,n;个数"<<endl;(n>m)中符号常量m"<<endl;(i=0;i<n;i++)"<<i+1<<"行:";"<<i+1<<"行:";equ(forcout<<x[i]<<"cout<<endl;return(0);}A,b, x,n);(i=0;i<n;i++)//x[n-2]=(Aug[n-2][n]-Aug[n-2][n-l]*x[nT])/Aug[n-2][n-2];for{s=0;fors=s+Aug[i][j]*x[j];cout<<"s="<<s<<endl;x[i]=(Aug[i][n]-s)/Aug[i][i];cout<<"x["<<i<<"]="<<x[i]<<endl;}2.3、高斯列主元程序调试对所编写的高斯列主元程序进行编译和链接，然后执行得如下所示的窗口，我们按命令输入增广矩阵的行数为3,输入3行5列的增广矩阵，运行界面为：丁面番辂人未轴數的不如T溝按顺序输入堆广矩阵斗12342.3、高斯列主元程序调试对所编写的高斯列主元程序进行编译和链接，然后执行得如下所示的窗口，我们按命令输入增广矩阵的行数为3,输入3行5列的增广矩阵，运行界面为：丁面番辂人未轴數的不如T溝按顺序输入堆广矩阵斗123451@4 6 92-0.11S38S3.9230&限方程组的解为=^[0J=-9■号9T9Q »LL1=58xE21=-34Pfc*eesainykeytocontinue—图2-2高斯列主元程序调试3、牛顿法解非线性方程组3.1、算法说明设Pk已知。第1步：计算函数F(P)二kf(p,q)f1(pk,qk)2kk-(3-1)第2步：计算雅可比矩阵—f(p,q)—f(p,q)j(p)_Qx i k k ay i k kk—f(p,q)—f(p,q)QX 2 k k ay 2 k k(3-2)第3步：求线性方程组J(P)AP=-F(P)k k第4步：计算下一点重复上述过程。P=重复上述过程。P=P+APk+1k(3-3)3.2、牛顿法解非线性方程组算法程序代码#include<iostream>#include<iomanip>#include<cmath>#defineN2//非线性方程组中方程个数、未知量个数#defineEpsilon0.0001//差向量1范数的上限#defineMax100//最大迭代次数usingnamespacestd;constintN2=2*N;intmain(){voidff(floatxx[N],floatyy[N]);voidffjacobian(floatxx[N],floatyy[N][N]);voidinv_jacobian(floatyy[N][N],floatinv[N][N]);voidnewdundiedai(floatx0[N],floatinv[N][N],floaty0[N],floatx1[N]);floatx0[N]={2.0,0.25},y0[N],jacobian[N][N],invjacobian[N][N],x1[N],errornorm;inti,j,iter=0;//for(i=0;i<N;i++)//cin>>x0[i];cout<<"初始解向量是："<<endl;for(i=0;i<N;i++)cout<<x0[i]<<"";cout<<endl;cout<<endl;do{iter=iter+1;cout<<"第"<<iter<<"次迭代开始"<<endl<<endl;//计算向量函数值---因变量向量存放在一维数组y0中，y0做函数ff的实参，ff(xO,yO)的调用将传出因变量向量ff(x0,y0);//计算雅克比矩阵，二维数组jacobian作为ffjacobian的实参;调用ffjacobian(x0,jacobian)，//按照地址传送方式，二维数组jacobian将传出雅克比矩阵的计算结果ffjacobian(x0,jacobian);//调用函数inv_jacobian计算雅克比矩阵的逆矩阵，计算结果存放在实参invjacobian中inv_jacobian(jacobian,invjacobian);//由近似解向量xO迭代计算近似解向量x1cout<〈〃第〃<<iter<<〃次迭代得到〃；newdundiedai(xO,invjacobian,yO,x1);//计算x1与xO的差向量的1范数errornorm=O；for(i=O；i<N；i++)errornorm=errornorm+fabs(x1[i]-xO[i])；if(errornorm<Epsilon)break；//为下次循环计算准备新的xOfor(i=O；i<N；i++)xO[i]=x1[i]；}while(iter<Max)；cin>>iter；returnO；}voidff(floatxx[N],floatyy[N]){floatx,y；inti；x=xx[O]；y=xx[1]；yy[O]=x*x-2*x+y+1；yy[1]=x*x-2*y*y+4；cout<<〃当前自变量向量为：〃<<endl;for(i=O；i<N；i++)cout<<setw(1O)<<setiosflags(ios::right)<<setprecision(5)<<setiosflags(ios::fixed)<<xx[i]<<〃 〃;cout<<endl;cout<<endl;cout<<〃当前因变量向量为：〃<<endl;for(i=O;i<N;i++)cout<<setw(10)<<setiosflags(ios::right)<<setprecision(5)<<setiosflags(ios::fixed)<<yy[i]<<"";cout<<endl;cout<<endl;}voidffjacobian(floatxx[N],floatyy[N][N]){floatx,y;inti,j;x=xx[0];y=xx[1];//jacobianhaven*nelementyy[0][0]=2*x-2;yy[0][1]=1;yy[1][0]=2*x;yy[1][1]=-4*y;cout<<"当前的雅克比矩阵是："<<endl;for(i=0;i<N;i++){for(j=0;j<N;j++)cout<<setw(10)<<setiosflags(ios::right)<<setprecision(5)<<setiosflags(ios::fixed)<<yy[i][j]<<" ";cout<<endl;}cout<<endl;}voidinv_jacobian(floatyy[N][N],floatinv[N][N]){floataug[N][N2],L;inti,j,k;cout<<"计算雅克比矩阵的逆矩阵"<<endl;for(i=0;i<N;i++){for(j=0;j<N;j++)aug[i][j]=yy[i][j];for(j=N;j<N2;j++)if(j==i+N)aug[i][j]=1;elseaug[i][j]=0;}for(i=0;i<N;i++){for(j=0;j<N2;j++)cout<<setw(10)<<setiosflags(ios::right)<<setprecision(5)<<setiosflags(ios::fixed)<<aug[i][j]<<"";cout<<endl;}cout<<endl;for(i=0;i<N;i++){for(k=i+1;k<N;k++){L=-aug[k][i]/aug[i][i];for(j=i;j<N2;j++)aug[k][j]=aug[k][j]+L*aug[i][j];}}for(i=0;i<N;i++){for(j=0;j<N2;j++)cout<<setw(10)<<setiosflags(ios::right)<<setprecision(5)<<setiosflags(ios::fixed)<<aug[i][j]<<"";cout<<endl;}cout<<endl;for(i=N-1;i>0;i--){for(k=i-1;k>=0;k--){L=-aug[k][i]/aug[i][i];for(j=N2-1;j>=0;j--)aug[k][j]=aug[k][j]+L*aug[i][j];}}for(i=0;i<N;i++){for(j=0;j<N2;j++)cout<<setw(10)<<setiosflags(ios::right)<<setprecision(5)<<setiosflags(ios::fixed)<<aug[i][j]<<"";cout<<endl;}cout<<endl;for(i=N-1;i>=0;i--)for(j=N2-1;j>=0;j--)aug[i][j]=aug[i][j]/aug[i][i];for(i=0;i<N;i++){for(j=0;j<N2;j++)cout<<setw(10)<<setiosflags(ios::right)<<setprecision(5)<<setiosflags(ios::fixed)<<aug[i][j]<<"";cout<<endl;for(j=N;j<N2;j++)inv[i][j-N]=aug[i][j];}cout<<endl;cout<〈"雅克比矩阵的逆矩阵是(InverseofJacobianis)"<<endl;for(i=0;i<N;i++){ for(j=0;j<N;j++)cout<<setw(10)<<setiosflags(ios::right)<<setprecision(5)<<setiosflags(ios::fixed)<<inv[i][j]<<" ";cout<<endl;}cout<<endl;}voidnewdundiedai(floatx0[N],floatinv[N][N],floaty0[N],floatx1[N]){inti,j;floatsum=0;for(i=0;i<N;i++){sum=0;for(j=0;j<N;j++)sum=sum+inv[i][j]*y0[j];龙贝格求积分算法xl[i]=xO[i]-sum;}cout<<"近似解向量"<<endl;for(i=0;i<N;i++)cout<<x1[i]<<" ";cout<<endl;cout<<endl;}#include<iostream>#include<cmath>#defineN2 //非线性方程组中方程个数、未知量个数#defineepsilon0.0001 //差向量1范数的上限#definemax10 //最大迭代次数usingnamespacestd;constintN2=2*N;intmain(){voidff(floatxx[N],floatyy[N]);voidffjacobian(floatxx[N],floatyy[N][N]);voidinv_jacobian(floatyy[N][N],floatinv[N][N]);voidnewdim(floatx0[N],floatinv[N][N],floaty0[N],floatx1[N]);floatx0[N]={2.0,0.25},y0[N],jacobian[N][N],invjacobian[N][N],x1[N],errornorm;inti,iter=0;cout<<"初始解向量："<<endl;for(i=0;i<N;i++)cout<<x0[i]<<" ";cout<<endl;do{iter=iter+1;cout<〈"第"<<iter<<"次迭代开始："<<endl;//jisuanxianglianghanshuzhi---yinbianliangxiangliangy0ff(x0,y0);//jisuanjacobianjuzhenjacobianffjacobian(x0,jacobian);//jisuanjacobianjuzhendenijuzheninvjacobianinv_jacobian(jacobian,invjacobian);//youjiexiangliangx0jisuanjiexiangliangx1newdim(x0,invjacobian,y0,x1);//jisuanchaxiangliangde1fanshuerrornorm=0;for(i=0;i<N;i++)errornorm=errornorm+fabs(x1[i]-x0[i]);if(errornorm<epsilon)break;for(i=0;i<N;i++)x0[i]=x1[i];}while(iter<max);return0;}voidff(floatxx[N],floatyy[N]){floatx,y;inti;x=xx[0];y=xx[1];//非线性方程组yy[0]=x*x-2*x-y+0.5;yy[1]=x*x+4*y*y-4;cout<〈"因变量向量："<<endl;for(i=0;i<N;i++)cout<<yy[i]<<"";cout<<endl;cout<<endl;}voidffjacobian(floatxx[N],floatyy[N][N]){floatx,y;inti,j;x=xx[0];y=xx[1];yy[0][0]=2*x-2;yy[0][1]=-1;yy[1][0]=2*x;yy[1][1]=8*y;cout<〈"雅克比矩阵："<<endl;for(i=0;i<N;i++){for(j=0;j<N;j++)cout<<yy[i][j]<<" ";cout<<endl;}cout<<endl;}voidinv_jacobian(floatyy[N][N],floatinv[N][N]){floataug[N][N2],L;inti,j,k;龙贝格求积分算法cout<〈"计算雅克比矩阵的逆："<<endl;for(i=0;i<N;i++){for(j=0;j<N;j++)aug[i][j]=yy[i][j];for(j=N;j<N2;j++)if(j==i+N)aug[i][j]=1;elseaug[i][j]=0;}for(i=0;i<N;i++){for(j=0;j<N2;j++)cout<<aug[i][j]<<"";cout<<endl;}cout<<endl;for(i=0;i<N;i++){for(k=i+1;k<N;k++){L=-aug[k][i]/aug[i][i];for(j=i;j<N2;j++)aug[k][j]=aug[k][j]+L*aug[i][j];}}for(i=0;i<N;i++){for(j=0;j<N2;j++)cout<<aug[i][j]<<"";cout<<endl;}cout<<endl;for(i=N-1;i>0;i--){for(k=i-1;k>=0;k--){L=-aug[k][i]/aug[i][i];for(j=N2-1;j>=0;j--)aug[k][j]=aug[k][j]+L*aug[i][j];}}for(i=0;i<N;i++){for(j=0;j<N2;j++)cout<<aug[i][j]<<"";cout<<endl;}cout<<endl;for(i=N-1;i>=0;i--)for(j=N2-1;j>=0;j--)aug[i][j]=aug[i][j]/aug[i][i];for(i=0;i<N;i++){for(j=0;j<N2;j++)cout<<aug[i][j]<<"";cout<<endl;for(j=N;j<N2;j++)inv[i][j-N]=aug[i][j];}cout<<endl;cout<<"雅克比矩阵的逆："<<endl;for(i=0;i<N;i++){for(j=0;j<N;j++)cout<<inv[i][j]<<"";cout<<endl;}cout<<endl;}voidnewdim(floatx0[N],floatinv[N][N],floaty0[N],floatx1[N]){inti,j;//计算非线性方程组的近似解向量floatsum=0;for(i=0;i<N;i++){sum=0;for(j=0;j<N;j++)sum=sum+inv[i][j]*y0[j];x1[i]=x0[i]-sum;}cout<<"近似解向量："<<endl;for(i=0;i<N;i++)cout<<x1[i]<<""<<endl;}3.3、牛顿法解非线性方程组算法程序调试

"C:\llsers\Administrator\De&ktop\Debug\Cpplexe'丈茹置向量是:2 0.25I叵II叵I2.00000 O.250OG当曲因变邑山量初：0.25000 0.25000当罪的雅兄比矩4是：2.00000-100000sers\Administrator\De&ktop\Debug\Cpplexe'sers\Administrator\De&ktop\Debug\Cpplexe'1.00000G.00QG00.0GS0G1.000000.25DQ0-0.500000.125000.25G00计算雅克比矩阵的逆矩阵2.00Q00-1.00Q001.D0GDGQ.DQGDQi+.GOSOO2.00S0G0.000001.000002.00Q00-1.000001.QOOQO0.Q00QO0.0000000000-2.000001.000002.O0QOO0.000000.50000O.250QO0.00000^.00000-2.000001.00000雅克比矩阵的逆矩阵是CInuerseofJacobianis)0.25000 0.12500-O.50QOO 0.25000第1叭钱厲.三刃诉乜：蛋向予1.30625 0.31250当前口变皇PI虽为］1.90625 0.31250当前凶变亘PI虽为：0.00S79 0.02441出前凶雅兒比距阵是：1.S1250 -1.000003.81250 2.50000计算雅克比矩阵的逆矩阵1.000000.90G90O.OQQOO1.000000.90G90O.OQQOO1.000003.31250 2.5O0QG应用本程序解方程组J一2x+y+1=°初始近似值xO,yO分别为2.00和0.25，经过100[x2-2y2+4=0次迭代求出X(l)=-1.58582和X(2)=0.58185。

4、龙贝格求积分算法4.1、算法说明生成J>K的逼近表R(J,K)，并以R(J+1,J+1)为最终解来逼近积分Ibf(x)dx沁R(J,J) (4—1)a逼近R(J,K)存在于一个特别的下三角矩阵中，第0列元素R(J,0)用基于2J个[a,b]子区间的连续梯形方法计算，然后利用龙贝格公式计算R(J,K)。当1<K<J时，第J行的兀素为R(J,K)=R(J,K-1)+R(J，K-1)-R(J-1，K-1) (4—2)4K—1当丨R(J,J)-R(J+1,J+l)l<tol时，程序在第(J+1)行结束。4.2、龙贝格求积分算法代码#include<iostream>#include<cmath>usingnamespacestd;#definef(x)sin(x*x)//举例函数#defineepsilon0.0001//精度#defineMAXREPT10//迭代次数，到最后仍达不到精度要求，则输出T(m=10).doubleRomberg(doubleaa,doublebb){//aa，bb积分上下限intm,n;//m控制迭代次数，而n控制复化梯形积分的分点数.n=2"mdoubleh，x;doubles，q;doubleep;//精度要求double*y二newdouble[MAXREPT];//为节省空间，只需一维数组//每次循环依次存储Romberg计算表的每行元素，以供计算下一行，算完后更新doublep;//p总是指示待计算元素的前一个元素(同一行)//迭代初值h=bb-aa;y[0]=h*(f(aa)+f(bb))/2.0;m=1;n=1;ep=epsilon+1.0;//迭代计算while((ep>=epsilon)&&(m<MAXREPT)){//复化积分公式求T2n(Romberg计算表中的第一列)，n初始为1,以后倍增p=0.0;for(inti=0;i<n;i++)//求Hnx=aa+(i+0.5)*h;p=p+f(x);}p=(y[0]+h*p)/2.0;//求T2n=l/2(Tn+Hn)，用p指示//求第m行元素，根据Romberg计算表本行的前一个元素(p指示),//和上一行左上角元素(y[k-1]指示)求得.s=1.0;for(intk=1;k<=m;k++){s=4.0*s;q=(s*p-y[k-1])/(s-1.0);y[k-1]=p;p=q;}p=fabs(q-y[m-1]);m=m+1;y[m-1]=q;n=n+n;h=h/2.0;}return(q);}intmain(){doublea，b;cout<〈"Romberg积分，请输入积分范围a,b:"<<endl;cin>>a>>b;cout<<"积分结果："<<Romberg(a,b)<<endl;return0;}4.3、龙贝格求积分算法程序调试我们以求解积分f(x)=sin(x2)，精度为0.0001，最高迭代10次为例，对所编写的龙贝格求积分算法程序进行编译和链接，经执行后得如下所示的窗口图4-2龙贝格求积分算法程序调试说明：应用Romberg算法求f(x)=sin(x2)在12〕区间上的精度为0.0001的积分为0.494508。图4-1算法流程图5、三次样条插值算法5.1三次样条插值算法说明表表5-1三次样条插值算法说明表策略描述包含m和m的方程0 N⑴三次紧压样条，确定S，（x），S《X）（如0 n果导数已知，这是“最佳选择”m=~(d-S'(x))-件0h0 0 20m=丄竹'(x)-d)-N h N N-1 2N-1（ii）natural三次样条（一条“松弛曲线”m=0，m=00 N（iii）外挂S"（x）到端点h(m—m)m=m—0 2 i0 i h1h (m -m)m=m—n-i n-i n-2N N-1 hN-2（iv）S"（x）是靠近端点的常量m=m m=m0 1' N N-1（v）在每个端点处指定S"（x）m=S''(x)m=S'(x)0 0 ' N N5.2、三次样条插值算法（压紧样条）程序调试我们将所编写的程序三次样条插值算法（压紧样条）程序进行调试勺二SSHY二Is:□勺二SSHY二Is:□耳琵4CBU曰俞1.X11-X1--ippP尸pF-P±n鹭❻二田±n£曰匚0±n>£A=iinWX±G=04±n =2±n址;2=2-OXTPU.-CzLnHo□o■■:一-D0■■:一-D0XPU.'t —寺晰兀i力一W-出三疾+¥条播值函数coLX±□吞=1-0_E立日回183百x邑.伍生陌咅回I些齐>Cka□22 O_<?QA^-^0!^+<：..：±>^3 *e_9®±^^ia1> — ® *tCs<0>2>^3 — —Xa3?G22*<>c* 3.33iS22->e<»: —=L>C23JS3-i- ——C23JS3-i- ——3= —图5-2三次样条插值算法程序运行结果（a）图5-2三次样条插值算法程序运行结果（b）3>-^3 ■+ 1_ —a3> —El,5^=131■_J S>图5-1三次样条插值算法（压紧样条）程序输入界面、运行结果色否堆续试验7：n亍匸口耳 hujjtocontiniuiu

运行结果分析:(5—1)0.06(x-1)+0.42(x-0)-0.06(x-1)+0.08(x-0)f(x)=J-0.42(2-x》-0.62(x-1)3-0.08(x-2)+2.62(x-1)(5—1)0.62(x-3)3+0.06(x-2)3-2.62(x-3)+1.44(x-2)作图程序(Matlab)：x1=0:0.01:1;y1=0.06*(x1—1).八3+0.42*(x1—0)「3—0.06*(x1—1)+0.08*(x1—0);x2=1:0.01:2;y2=—0.42*(x2—2)「3—0.62*(x2—1).八3—0.08*(x2—2)+2.62*(x2—1);x3=2:0.01:3;y3=0.62*(x3—3).八3+0.06*(x3—2)「3—2.62*(x3—3)+1.44*(x3—2);X=[0123];Y=[00.521.5];plot(x1,y1,x2,y2,x3,y3,X,Y,'*')gtext('S1')gtext('S2')gtext('S3')图形为:图形为:图5-3三次样条插值算法Matlab作图分析5.3、三次样条插值算法(压紧样条)代码#include<iostream>#include<iomanip>usingnamespacestd;constintmax=50;floatx[max],y[max],h[max];floatc[max],a[max],fxym[max];floatf(intx1,intx2,intx3){floata=(y[x3]-y[x2])/(x[x3]-x[x2]);floatb=(y[x2]-y[x1])/(x[x2]-x[x1]);return(a-b)/(x[x3]-x[x1]);}//求差分voidcal_m(intn){ //用追赶法求解出弯矩向量M……floatB[max];B[0]=c[0]/2;for(inti=1;i<n;i++)B[i]=c[i]/(2-a[i]*B[i-1]);fxym[0]=fxym[0]/2;for(i=1;i<=n;i++)fxym[i]=(fxym[i]-a[i]*fxym[i-1])/(2-a[i]*B[i-1]);for(i=n-1;i>=0;i--)fxym[i]=fxym[i]-B[i]*fxym[i+1];

voidprintout(intn);intmain(){intn,i;charch;do{cout<<"输入x的最大下标:"；cin>>n;for(i=0；i<=n；i++){cout<<"PleaseputinX"<<i<<':'；cin>>x[i]；cout<<"PleaseputinY"<<i<<':'；cin>>y[i]；}for(i=0；i<n；i++) //求步长h[i]=x[i+1]-x[i]；:两端的cout<〈"Please输入边界条件\n1 :已知两端的一阶导数\n2:两端的二阶导数已知\n默认：自然边界条件\n";intt；floatf0,f1;cin>>t;switch(t){case1:cout<<"输入Y0\'Y"<<n<<"\'\n";cin>>f0>>f1;c[0]=1;a[n]=1;fxym[0]=6*((y[1]-y[0])/(x[1]-x[0])-f0)/h[0];fxym[n]=6*(f1-(y[n]-y[n-1])/(x[n]-x[n-1]))/h[n-1];break;case2:cout<<"输入Y0\"Y"<<n<<"\"\n";cin>>f0>>f1;c[0]=a[n]=0;fxym[0]=2*f0;fxym[n]=2*f1;break;default:cout<<"不可用\n";//待定};for(i=1;i<n;i++)fxym[i]=6*f(i-1,i,i+1);for(i=1;i<n;i++){a[i]=h[i-1]/(h[i]+h[i-1]);c[i]=1-a[i];}a[n]=h[n-1]/(h[n-1]+h[n]);cal_m(n);cout<<"\n输出三次样条插值函数：\n";printout(n);cout<<"Doyouhaveanothertry?y/n:";cin>>ch;}while(ch=='y' ||ch=='Y');return0;}voidprintout(intn){cout<<setprecision(6);for(inti=0;i<n;i++){cout<<i+1<<":["<<x[i]<<","<<x[i+1]<<"]\n"<<"\t";cout<<"S"<<i+1<<"=";floatt=fxym[i]/(6*h[i]);if(t>0)cout<〈—1<<"*(x—"<<x[i+l]<〈"厂3";elsecout<<-t<<"*(x—"<<x[i+l]<〈"厂3";t=fxym[i+1]/(6*h[i]);if(t>0)cout<<"+"<<t<<"*(x—"<<x[i]<〈"厂3";elsecout<<"-"<<t<<"*(x—"<<x[i]<〈"厂3";cout<<"\n\t";t=(y[i]—fxym[i]*h[i]*h[i]/6)/h[i];if(t>0)cout<<"—"<<t<<"*(x—"<<x[i+1]<<")";elsecout<<"—"<<—t<<"*(x—"<<x[i+1]<<")";t=(y[i+1]—fxym[i+1]*h[i]*h[i]/6)/h[i];if(t>0)cout<<"+"<<t<<"*(x—"<<x[i]<<")";elsecout<<"—"<<—t<<"*(x—"<<x[i]<<")";cout<<endl;}cout<<endl;}

6、M次多项式曲线拟合6.1、算法说明设“Xk'yk)l有N个点，横坐标是确定的。最小二乘抛物线的系数表示为y=f(x)=Ax2+y=f(x)=Ax2+Bx+C求解A，B和C的线性方程组为A+(6—1)乙X4k、k=l 丿乙X3k、k=l 丿乙X2k、k=l 丿(y )乙X3k、k=1 丿乙X2k、k=1 丿送X2]c=^EyJkB+送x」C=^yIk=1k丿'送x]B+NC=^ykIk=1丿X2kkk=1Xk丿 kk'k=1 丿 k=1kk=1(6—2)6.2、M次多项式曲线拟合算法流程图幵胎P利用吕斯列主元法计芹A^=F'*FaB-Fr^VrM-*i<—i-l-1图6-1算法流程图6.3、M次多项式曲线拟合算法程序调试6.3、我们按命令依次输入命令如下命令后，得程序执行结果如下

驍际沁的个如请输人4个点的X坐标序列:0123请输入4个点的¥坐标序列:00.521_5请输入需要拟合的次如拟合后的锹多项式系数划幕次.由高到低：0 0 -0_9 -0_25 1.35 -0.15拟合后的锹多项式为二P<«>=+0*xA5+0*x^+-0・9*x^3+-B・25 *1・35 *-0・15Pi-essani/keytocontInue图6-2算法调试图作图程序：X=[-3-113];Y=[15515];x=-3:0.01:3;y=0.875*x「2-1.7*x+2.145;plot(X,Y,'go',x,y,'r-')gtext('拟合曲线')6.4、 M次多项式曲线拟合算法代码#include<iostream>#include<cmath>#defineMAX20usingnamespacestd;图形为:图6-3图形分析//求解任意可逆矩阵的逆，X为待求解矩阵，E图形为:图6-3图形分析voidinv(doubleX[MAX][MAX],intn,doubleE[MAX][MAX]){inti,j,k;doubletemp=0;for(i=0;i<MAX;i++){for(j=0;j<MAX;j++)if(i==j)E[i][j]=1;}for(i=0;i<n-1;i++){temp=X[i][i];for(j=0;j<n;j++){X[i][j]=X[i][j]/temp;E[i][j]=E[i][j]/temp;}for(k=0;k<n;k++){if(k==i)continue;temp=-X[i][i]*X[k][i];for(j=0;j<n;j++)X[k][j]=X[k][j]+temp*X[i][j];E[k][j]=E[k][j]+temp*E[i][j];}}}}intmain(){intn,M,i,j,k;doubleX[MAX]={0},Y[MAX]={0},F[MAX][MAX]={0},B[MAX]={0};doubleA[MAX][MAX]={0},BF[MAX][MAX]={0},C[MAX]={0},E[MAX][MAX]={0};cout<<"M次多项式曲线拟合\n请先输入待拟合的点的个数:"；cin>>n;cout<<"\n请输入"<<n<〈"个点的X坐标序列:\n";for(i=0；i<n；i++)cin>>X[i];cout<<"\n请输入"<<n<〈"个点的Y坐标序列:\n";for(i=0;i<n;i++)cin>>Y[i];cout<<"\n请输入需要拟合的次数:";cin>>M;for(i=0;i<n;i++)for(k=1;k<=M+1;k++)F[i][k-1]=pow(X[i],k-1);//求F的转置for(i=0;i<n;i++){for(j=0;j<M+1;j++){BF[j][i]=F[i][j];}}//计算其转置的BF与F的乘for(i=0;i<M+1;i++)for(j=0;j<M+1;j++)for(k=0;k<n;k++)A[i][j]+=BF[i][k]*F[k][j];//计算F的转置BF与Y的乘for(i=0;i<M+1;i++)for(j=0;j<n;j++)B[i]+=BF[i][j]*Y[j];//调用inv函数求解矩阵A的逆矩阵Einv(A,n,E);//计算A的逆BF与B的乘for(i=0;i<M+1;i++)for(j=0;j<n;j++)C[i]+=E[i][j]*B[j];cout<<"\n拟合后的"<<M<〈"次多项式系数为，幂次由高到低：\n";for(i=M;i>=0;i--)cout<<C[i]<<"\t";cout<<"\n拟合后的"<<M<〈"次多项式为:\n";cout<<"\nP(x)=";for(i=M;i>=0;i--){if(i==0)cout<<"+"<<C[i];elsecout<<〃+〃<<C[i]<<〃*x""<<i;}cout<<endl;return0;}7、牛顿-拉弗森迭代解非线性方程7.1、牛顿-拉弗森迭代解非线性方程算法概要使用初始近似值P，利用迭代P=P-f卩丘i)其中k=1,2,...计算函数f(x)=0的根的0 k k-1f(p) 八k-1近似值。7.2、牛顿-拉弗森迭代解非线性方程算法程序调试图7-1牛顿-拉弗森迭代解非线性方程算法程序调试7.3、牛顿-拉弗森迭代解非线性方程算法代码#include<iostream>usingnamespacestd;#include<cmath>doublef(doublex){returnx*x—2*xT}doubledf(doublex){return2*x—2;}intmain(){intk=0,max1=0;doublep0,delta,epsilon,p1,err,relerr,y;cout<〈"牛顿拉弗森法解非线性方程f(x)=x^2—2x—1"<<endl;cout<<"初始值为p0=4.8"<<endl;p0=4.&delta=0.0001,epsilon=0.0001,max1=100;for(k=1;k<=max1;k++){p1=p0—f(p0)/df(p0);err二fabs(p1-p0);relerr=2*err/(fabs(p1)+delta);p0=p1;牛顿-拉弗森迭代解非线性方程if((err<delta)||(relerr<delta)||(fabs(y)<epsilon))break;}cout<〈"方程的根的近似值为:"<<pO<<endl;cout<〈"方程的根的误差估计为:"<<err<<endl;cout<<"迭代次数为:"<<k<<endl;cout<〈"方程在p0点的函数值f(pO):"<<y<<endl;return0;}8、不动点法解非线性方程8.1、 算法说明先将f（x）改写成x=g（x）然后对x=g（x）进行迭代，即X=g（xk-1）其中k=12…然后判断1Xk-Xk-il<£是否成立，成立则返回Xk，不成立就重复以上步骤8.2、 不动点法解非线性方程算法程序调试我们将编写好的不动点迭代法解非线性方程算法程序进行编译，链接和执行后得如下所示结果。|不动点法解非线性方程F3W2忽片稚社［①，1】上有齡，初始值>^Jp0=0用动点的近似值知2憧序选代次数为活近似值的误差为=G-10342e-005障解不动点近似值的序列二0 1 1-75 1.98438 1.99994PressanykuytocontinueH图8-1程序调试图8.3、不动点法解非线性方程算法代码#include<iostream>#include<cmath>#include<iomanip>#defineMAX20#defineepsle—10usingnamespacestd;doubleg(doublex){return1—x*x/4+x;}main(){doubleP[MAX]二{0},err=0.0,relerr=0.0,tol=0.0,p=0.0,p0=0.0;intk=0,max1=0,i=0;cout<〈"不动点法解非线性方程f(x)=l-x"2/2"<<endl;cout<<〃方程在[0,1]上有解，初始值为p0=0〃<<endl;//初始化P[0]=p0=0;max1=100;tol=0.001;for(k=2;k<=max1;k++){P[k-1]=g(P[k-2]);err=fabs(P[k-1]-P[k-2]);relerr=err/(fabs(P[k-1]+eps));p=P[k-1];if((err<tol)||(relerr<tol))break;}if(k==max1)cout<<〃迭代次数超过允许的最大迭代次数！"<<endl;cout<<〃不动点的近似值为:〃<<p<<endl;cout<<〃程序迭代次数为:"<<k<<endl;cout<〈〃近似值的误差为:"<<err<<endl;cout<<〃求解不动点近似值的序列:"<<endl;for(i=0;i<k;i++){ cout<<setw(10)<<P[i];}cout<<endl;return0;9.拉格朗日插值9.1、算法说明有N+1个点(x,y),(x,y (x,y)的次数最高为N的多项式P(x)的构造方法，它具有0 0 11 nn NP(x)=£yL (x)的形式，其中L (x)是基于节点N kN,k N,kk=0(xL (x)=l-x)(x-x)・・・(x-x )(x-x)・・・(x-x)气( 3丫(n\的拉格朗日系数多项式x-xx-xx-x x-x丿k1 k k-1 k k+1 knN，k (X-Xk0对每个固定的k,拉格朗日系数多项式L (x)具有性质为:N,k9.2、拉格朗日插值算法程序调试1,k=i0,k丰i(9—1)首先编写好程序，然后编译链接，运行程序，按照程序提示依次输入点的个数、点数对应的x值、y值、纵坐标、再输入x值，输入结果如下：癩点的个数门柑请输入=3请输4请输入曲須盂请输入昭 图9-1拉格朗日插值算法程序调试厦过拉格朋日插值公式丙得的结果：当兀勻时，尸§Pressanykeytocontinu.e图9-2拉格朗日插值算法程序计算结果9.3、 拉格朗日插值算法程序代码：#include〈iostream>usingnamespacestd;doublefunc(doubleX,intk,doublex[],intn);intmain(){doubleSn=0;intn;cout〈〈"请输入点的个数n:";cin>>n;double*x=(double*)malloc(n*sizeof(double));double*y=(double*)malloc(n*sizeof(double));doubleX;inti;for(i=0;i<n;i++){cout<<"请输入x"<<i+l<<",y"<<i+l<〈"："<<endl;cin>>x[i]>>y[i];}cout<<"请输入x";cin>>X;for(i=0;i<