经典四阶龙格库塔法解一阶微分方程组(共45页)

1、精选优质文档-倾情为你奉上1.经典四阶龙格库塔法解一阶微分方程组1.1运用四阶龙格库塔法解一阶微分方程组算法分析（1-1）， （1-2） （1-3） （1-4）（1-5）（1-6）（1-7）（1-8）（1-9） （1-10）经过循环计算由 推得 每个龙格-库塔方法都是由一个合适的泰勒方法推导而来，使得其最终全局误差为,一种折中方法是每次进行若干次函数求值，从而省去高阶导数计算。4阶龙格-库塔方法(RK4)是最常用的，它适用于一般的应用，因为它非常精准，稳定，且易于编程。1.2经典四阶龙格库塔法解一阶微分方程流程图 图1-1 经典四阶龙格库塔法解一阶微分方程流程图1.3经典四阶龙格库塔法解一阶微

2、分方程程序代码： #include <iostream>#include <iomanip>using namespace std;void RK4( double (*f)(double t,double x, double y),double (*g)(double t,double x, double y) ,double initial3, double resu3,double h) double f1,f2,f3,f4,g1,g2,g3,g4,t0,x0,y0,x1,y1; t0=initial0;x0=initial1;y0=initial2; f1=f(

3、t0,x0,y0); g1=g(t0,x0,y0); f2=f(t0+h/2, x0+h*f1/2,y0+h*g1/2); g2=g(t0+h/2, x0+h*f1/2,y0+h*g1/2); f3=f(t0+h/2, x0+h*f2/2,y0+h*g2/2); g3=g(t0+h/2, x0+h*f2/2,y0+h*g2/2); f4=f(t0+h, x0+h*f3,y0+h*g3); g4=g(t0+h, x0+h*f3,y0+h*g3); x1=x0+h*(f1+2*f2+2*f3+f4)/6; y1=y0+h*(g1+2*g2+2*g3+g4)/6;resu0=t0+h;resu1=x

4、1;resu2=y1;int main() double f(double t,double x, double y); double g(double t,double x, double y); double initial3,resu3; double a,b,H; double t,step; int i; cout<<"输入所求微分方程组的初值t0,x0,y0:" cin>>initial0>>initial1>>initial2; cout<<"输入所求微分方程组的微分区间a,b:"

5、 cin>>a>>b; cout<<"输入所求微分方程组所分解子区间的个数step:" cin>>step; cout<<setiosflags(ios:right)<<setiosflags(ios:fixed)<<setprecision(10); H=(b-a)/step; cout<< initial0<<setw(18)<<initial1<<setw(18)<<initial2<<endl; for(i=0;

6、i<step;i+) RK4( f,g ,initial, resu,H); cout<<resu0<<setw(20)<<resu1<<setw(20)<<resu2<<endl; initial0=resu0;initial1=resu1;initial2=resu2; return(0);double f(double t,double x, double y) double dx; dx=x+2*y;return(dx);double g(double t,double x, double y)double

7、dy;dy=3*x+2*y;return(dy);1.4经典四阶龙格库塔法解一阶微分方程程序调试结果图示： 应用所编写程序计算所给例题： 其中初值为求解区间为0,0.2。 计算结果为: 图1-2 经典四阶龙格库塔法解一阶微分方程算法程序调试图专心-专注-专业2.高斯列主元法解线性方程组2.1高斯列主元法解线性方程组算法分析使用伪代码编写高斯消元过程：for k=1 to n-1 dofor i=k+1 to nl<=a(i,k)/a(k,k)for j=k to n doa(i,j)<=a(i,j)-l*a(k,j) end %end of for jb(i)<=b(i)-l

8、*b(k) end %end of for iend %end of for k最后得到A，b可以构成上三角线性方程组接着使用回代法求解上三角线性方程组 因为高斯消元要求a（k,k）0（k=1,2,3n-1）这就需要对高斯消元过程进行完善，即使用高斯列主元法：其步骤为： 找主元：计算 ,并记录其所在行r , 交换第r行与第k行； 以第k行为工具行处理以下各行，使得从第k列的第k+1行到第n行的元素全部为0； 得到增广矩阵的上三角线性方程组；使用回代法对上三角线性方程组进行求解2.2高斯列主元法解线性方程组流程图FFTT开始输入A交换A中r、k两行输出x结束图2-1 高斯列主元法解线性方程组流程

9、图2.3高斯列主元法解线性方程组程序代码#include<iostream>#include <cmath>#define N 3using namespace std;void main()int i,j,k,n,p;float t,s,m,aNN,bN,xN;cout<<"请输入方程组的系数"<<endl;for(i=0;i<N;i+)for(j=0;j<N;j+)cin>>aij;cout<<"请输入方程组右端的常数项:"<<endl;for(i=0;i

10、<N;i+) cin>>bi;for(j=0;j<N-1;j+) p=j; for(i=j+1;i<N;i+) /寻找系数矩阵每列的最大值 if(fabs(aij)>fabs(apj) p=i; if(p!=j) /交换第p行与第j行for(k=0;k<N;k+)t=ajk;ajk=apk;apk=t; /交换常数项的第p行与第j行 t=bp; bp=bj; bj=t; /把系数矩阵第j列ajj下面的元素变为0for(i=j+1;i<N;i+) m=-aij/ajj; for(n=j;n<N;n+) ain=ain+ajn*m;bi=bi+

11、bj*m; /求方程组的一个解xN-1=bN-1/aN-1N-1; /回代法求方程组其他解for(i=N-2;i>=0;i-) s=0.0;for(j=N-1;j>i;j-)s=s+aij*xj; xi=(bi-s)/aii; cout<<"方程组的解如下:"<<endl;for(i=0;i<=N-1;i+) cout<<xi<<endl;2.4高斯列主元法解线性方程组程序调试结果图示：求解下列方程组 图2-2 高斯列主元法解线性方程组程序算法调试图 3.牛顿法解非线性方程组3.1牛顿法解非线性方程组算法说明

12、 牛顿法主要思想是用 进行迭代 。因此首先需要算出的雅可比矩阵，再求过求出它的逆，当它达到某个精度时即停止迭代。 具体算法如下：首先设已知。：计算函数 （3-1） ：计算雅可比矩阵 （3-2） （3-3）：求线性方程组的解。：计算下一点重复上述过程。3.2牛顿法解非线性方程组流程图图3-1 牛顿法解非线性方程组流程图3.3牛顿法解非线性方程组程序代码#include<iostream>#include<cmath>#define N 2 / 非线性方程组中方程个数、未知量个数 #define Epsilon 0.0001 / 差向量1范数的上限#define Max 1

13、00 /最大迭代次数using namespace std;const int N2=2*N;int main()void ff(float xxN,float yyN);/计算向量函数的因变量向量yyNvoid ffjacobian(float xxN,float yyNN);/计算雅克比矩阵yyNNvoid inv_jacobian(float yyNN,float invNN);/计算雅克比矩阵的逆矩阵invvoid newdundiedai(float x0N, float invNN,float y0N,float x1N);/由近似解向量 x0 计算近似解向量 x1float x0

14、N=2.0,0.25,y0N,jacobianNN,invjacobianNN,x1N,errornorm;int i,j,iter=0;/如果取消对x0的初始化，撤销下面两行的注释符，就可以由键盘向x0读入初始近似解向量/for( i=0;i<N;i+)/ cin>>x0i;cout<<"初始近似解向量："<<endl;for (i=0;i<N;i+)cout<<x0i<<" " cout<<endl;cout<<endl;doiter=iter+1;cou

15、t<<"第 "<<iter<<" 次迭代开始"<<endl;/计算向量函数的因变量向量 y0ff(x0,y0);/计算雅克比矩阵 jacobianffjacobian(x0,jacobian);/计算雅克比矩阵的逆矩阵 invjacobianinv_jacobian(jacobian,invjacobian);/由近似解向量 x0 计算近似解向量 x1newdundiedai(x0, invjacobian,y0,x1);/计算差向量的1范数errornormerrornorm=0;for (i=0;i&l

16、t;N;i+)errornorm=errornorm+fabs(x1i-x0i);if (errornorm<Epsilon) break;for (i=0;i<N;i+)x0i=x1i; while (iter<Max);return 0;void ff(float xxN,float yyN)float x,y; int i; x=xx0;y=xx1;yy0=x*x-2*x-y+0.5;yy1=x*x+4*y*y-4; cout<<"向量函数的因变量向量是： "<<endl;for( i=0;i<N;i+) cout<

17、;<yyi<<" " cout<<endl;cout<<endl;void ffjacobian(float xxN,float yyNN) float x,y; int i,j;x=xx0;y=xx1;/jacobian have n*n elementyy00=2*x-2;yy01=-1;yy10=2*x;yy11=8*y; cout<<"雅克比矩阵是： "<<endl;for( i=0;i<N;i+) for(j=0;j<N;j+) cout<<yyij<

18、;<" "cout<<endl; cout<<endl;void inv_jacobian(float yyNN,float invNN)float augNN2,L; int i,j,k; cout<<"开始计算雅克比矩阵的逆矩阵 ："<<endl; for (i=0;i<N;i+) for(j=0;j<N;j+) augij=yyij; for(j=N;j<N2;j+)if(j=i+N) augij=1;else augij=0; for (i=0;i<N;i+) for(

19、j=0;j<N2;j+) cout<<augij<<" " cout<<endl;cout<<endl;for (i=0;i<N;i+) for (k=i+1;k<N;k+) L=-augki/augii;for(j=i;j<N2;j+) augkj=augkj+L*augij; for (i=0;i<N;i+) for(j=0;j<N2;j+) cout<<augij<<" " cout<<endl;cout<<endl;

20、 for (i=N-1;i>0;i-) for (k=i-1;k>=0;k-) L=-augki/augii;for(j=N2-1;j>=0;j-) augkj=augkj+L*augij; for (i=0;i<N;i+) for(j=0;j<N2;j+) cout<<augij<<" " cout<<endl;cout<<endl;for (i=N-1;i>=0;i-) for(j=N2-1;j>=0;j-)augij=augij/augii;for (i=0;i<N;i+)

21、 for(j=0;j<N2;j+) cout<<augij<<" " cout<<endl; for(j=N;j<N2;j+) invij-N=augij;cout<<endl;cout<<"雅克比矩阵的逆矩阵： "<<endl;for (i=0;i<N;i+) for(j=0;j<N;j+) cout<<invij<<" " cout<<endl;cout<<endl;void newdun

22、diedai(float x0N, float invNN,float y0N,float x1N)int i,j;float sum=0;for(i=0;i<N;i+) sum=0; for(j=0;j<N;j+) sum=sum+invij*y0j;x1i=x0i-sum; cout<<"近似解向量："<<endl;for (i=0;i<N;i+) cout<<x1i<<" " cout<<endl;cout<<endl;3.4牛顿法解非线性方程组程序调试结果图

23、示：图3-2 牛顿法解非线性方程组程序算法调试图图3-3 牛顿法解非线性方程组程序算法调试图 图3-4 牛顿法解非线性方程组程序算法调试图4.龙贝格求积分算法4.1龙贝格求积分算法分析生成j>=k 的近似积分结果逼近表，并以R(j+1,j+1)为最终解来逼近积分。R(j,0)=T(j), j>=0,T(j)为区间逐次减半递推梯形求积分公式算出的结果；R(j,1)=S(j), j>=1,S(j)为区间逐次减半递推辛普森求积分公式算出的结果；（4-1）R(j,2)=B(j), j>=2,B(j)为递推布尔求积分公式算出的结果； （4-2）生成的逼近表，并以为最终解来逼近积分

24、（4-3）逼近存在于一个特别的下三角矩阵中，第0列元素用基于个a,b子区间的连续梯形方法计算，然后利用龙贝格公式计算。当时，第行的元素为当时，程序在第行结束。4.2龙贝格积分算法流程图 图4-1 龙贝格积分算法流程图4.3龙贝格积分算法程序代码#include<iostream>using namespace std;#include<stdio.h>#include<math.h>#define f(x) (sin(x) /列举函数#define N_H 20#define MAX 10 /最大迭代次数#define a 0 /所求积分的上下限#defin

25、e b 1#define epsilon 0.0001 /所需精度double Romberg(double aa,double bb,long int n) int i; double sum,h=(bb-aa)/n;sum=0; for(i=1;i<n;i+) sum+=f(aa+i*h); sum+=(f(aa)+f(bb)/2; return (h*sum);void main() int i; long int n=N_H,m=0; double T2MAX+1; T10=Romberg(a,b,n); n*=2; for(m=1;m<MAX;m+) for(i=0;i&

26、lt;m;i+) T0i=T1i; T10=Romberg(a,b,n); n=n*2; for (i=1;i<=m;i+) T1i=T1i-1+(T1i-1-T0i-1)/(pow(2,2*m)-1); if(T0m-1-T1m)<epsilon) cout<<"T="<<T1m<<endl;return; 4.4龙贝格积分算法调试图求在区间上的精度为0.0001的积分 图4-2 龙贝格积分程序算法调试图5.三次样条插值算法5.1三次样条插值基本算法说明： 表5-1三次样条插值基本算法说明策略描述包含和的方程(i)三次紧压样

27、条，确定，（如果导数已知，这是“最佳选择”）(ii)natural三次样条（一条“松弛曲线”），(iii)外挂到端点若已知N+1个点的及其一阶导数的边界条件S(a)=和S(b)=,则存在唯一的三次样条曲线。构造并求解下列线性方程组（5-1） （5-2）当得到系数 后，可利用如下公式计算样条函数 （5-3） 为了更有效的计算，每个三次多项式 可表示成嵌套乘的形式，也可以写为如下形式：（5-4）5.2三次样条插值算法程序代码#include<stdio.h> #include<stdlib.h> #define MAX 4 double *diff(double X,int

28、 n) int i=0; double *H=NULL; H=(double*)malloc(n-1)*sizeof(double); for(i=1;i<=n-1;i+) Hi-1=Xi-Xi-1; return H; double *divide(double Y,int N,double H) int i=0; double *D=NULL; D=(double*)malloc(N*sizeof(double); for(i=0;i<N;i+) Di=Yi/Hi; return D; main() doubleXMAX=0,1,2,3,YMAX=0,0.5,2.0,1.5,S

29、MAXMAX=0,temp=0.0,MMAX=0; int N=MAX-1,i=0,k=0; double AMAX-1-2=0,BMAX-1-1=0,CMAX-1-1=0; double dx0=0.2,dxn=1.0; double *H=NULL,*D=NULL,*U=NULL; printf("求解过点（0，0.0），（1，0.5），（2，2.0）和（3，1.5）n而且一阶导数边界条件S'(0)=0.2和S'(2)=-1的三次压紧样条曲线nn"); H=diff(X,MAX); D=divide(diff(Y,MAX),N,H); for(i=1;i

30、<N-2;i+) Ai=Hi+1; for(i=0;i<N-1;i+) Bi=2*(Hi+Hi+1); for(i=1;i<N-1;i+) Ci=Hi+1; U=diff(D,N); for(i=0;i<N;i+) Ui=Ui*6; B0=B0-H0/2; U0=U0-3*(D0-dx0); BN-2=BN-2-HN-1/2; UN-2=UN-2-3*(dxn-DN-1); for(k=2;k<=N-1;k+) temp=Ak-2/Bk-2; Bk-1=Bk-1-temp*Ck-2; Uk-1=Uk-1-temp*Uk-2; MN-1=UN-2/BN-2; for

31、(k=N-2;k>=1;k-) Mk=(Uk-1-Ck-1*Mk+1)/Bk-1; M0=3*(D0-dx0)/H0-M0/2; MN=3*(dxn-DN-1)/HN-1-MN-1/2; for(k=0;k<=N-1;k+) Sk0=(Mk+1-Mk)/(6*Hk); Sk1=Mk/2; Sk2=Dk-Hk*(2*Mk+Mk+1)/6; Sk3=Yk; printf("求得的三次紧压样条曲线的矩阵S为:n"); for(i=0;i<MAX-1;i+) for(k=0;k<MAX;k+) printf("%lft",Sik); pr

32、intf("n"); for(i=0;i<N;i+) printf(y=%+lfx3%+lfx2%+lfx%+lf",Si0,Si1,Si2,Si3); getchar(); 5.3三次样条插值算法调试运行结果求解三次紧压样条曲线，以过点（0，0.0），（1，0.5），（2，2.0）和（3，1.5），而且一阶导数边界条件S'(0)=0.2和S'(2)=-1的三次压紧样条曲线。图5-1 三次样条插值程序算法调试图5.4三次样条插值用matlab画图>> x1=0:.01:1; y1=polyval(S(1,:),x1-X(1);x2

33、=1:.01:2; y2=polyval(S(2,:),x2-X(2);x3=2:.01:3; y3=polyval(S(2,:),x3-X(3);plot(x1,y1,x2,y2,x3,y3,X,Y,'.')运行结果如下： 图5-2 三次样条插值作图 6.M次多项式曲线拟合6.1 M次多项式曲线拟合算法说明：设有N个点，横坐标是确定的。最小二乘法抛物线系数表示为，（6-1）求解A，B和C的线性方程组为（6-2）整体上考虑近似函数同所给数据点(i=0,1,m)误差 (i=0,1,m)的大小，常用的方法有以下三种：一是误差(i=0,1,m)绝对值的最大值，即误差向量的范数；二是误

34、差绝对值的和，即误差向量r的1范数；三是误差平方和的算术平方根，即误差向量r的2范数；前两种方法简单、自然，但不便于微分运算 ，后一种方法相当于考虑 2范数的平方，因此在曲线拟合中常采用误差平方和来 度量误差 (i=0，1，m)的整体大小。数据拟合的具体作法是：对给定数据 (i=0,1,，m)，在取定的函数类中,求,使误差(i=0,1,m)的平方和最小。6.2 M次多项式曲线拟合流程图FT开始输入x0、输出极小值点和极小值结束图6-1 M次曲线多项式曲线拟合流程图6.3 M次多项式曲线拟合算法程序代码#include<stdio.h> #include<math.h>

35、#define MAX 20 void inv(double XMAXMAX,int n,double EMAXMAX) int i=0,j=0,k=0; double temp=0.0; for(i=0;i<MAX;i+) for(j=0;j<MAX;j+) if(i=j) Eij=1; for(i=0;i<n-1;i+) temp=Xii; for(j=0;j<n;j+) Xij=Xij/temp; Eij=Eij/temp; for(k=0;k<n;k+) if(k=i) continue; temp=-Xii*Xki; for(j=0;j<n;j+)

36、 Xkj=Xkj+temp*Xij; Ekj=Ekj+temp*Eij; main() int n=0,M=0,i=0,j=0,k=0; double XMAX=0,YMAX=0,FMAXMAX=0,BMAX=0; double AMAXMAX=0,BFMAXMAX=0,EMAXMAX=0,CMAX=0; printf("tttM次多项式曲线拟合nn输入拟合的点的个数:"); scanf("%d",&n); printf("n输入%d个点的X坐标序列:n",n); for(i=0;i<n;i+) scanf("

37、%lf",&Xi); printf("n输入%d个点的Y坐标序列:n",n); for(i=0;i<n;i+) scanf("%lf",&Yi); printf("n输入拟合次数:"); scanf("%d",&M); for(i=0;i<n;i+) for(k=1;k<=M+1;k+) Fik-1=pow(Xi,k-1); for(i=0;i<n;i+) for(j=0;j<M+1;j+) BFji=Fij; for(i=0;i<M+1;i+)

38、 for(j=0;j<M+1;j+) for(k=0;k<n;k+) Aij+=BFik*Fkj; for(i=0;i<M+1;i+) for(j=0;j<n;j+) Bi+=BFij*Yj; inv(A,n,E); for(i=0;i<M+1;i+) for(j=0;j<n;j+) Ci+=Eij*Bj; printf("n拟合后的%d次多项式系数：n",M); for(i=M;i>=0;i-) printf("%lft",Ci); printf("n拟合后的%d次多项式为:n",M); p

39、rintf("nP(x)="); for(i=M;i>=0;i-) if(i=0) printf("%+.3lfn",Ci); else printf("%+.3lf*x%d",Ci,i); getchar(); 6.4 M次多项式曲线拟合算法程序调试结果：求解（-3，3），（0，1），（2，1），（4，3）四个点拟合2次后的2后的多项式图6-2 M次曲线多项式曲线拟合程序算法调试图6.5 M次多项式曲线拟合x=-3:0:2:4;y=3:1:1:3;plot(x,y,'b')图6-3 M次多项式曲线拟合画图7.二

40、分法解非线性方程7.1 二分法解非线性方程算法说明：.是起始区间，是中点。.是第二个区间，它包含零点r，同时是中点，区间的宽度范围是的一半。.得到第n个区间（包含r，并有中点）后，可构造出，它也包括r，宽度范围是的一半。 7.2二分法解非线性方程程序代码#include<iostream>using namespace std;#define eps 0.001double f(double x) return 2-x*x;int main()double ga,gb,gc,a,b,c; cout<<" *"<<endl;cout<

41、<" 用二分法寻找方程2-x*x=0的根"<<endl;cout<<" 求根区间为(1,2)"<<endl; cout<<" *"<<endl;a=1,b=2;ga=f(a);gb=f(b);while(b-a)>eps)c=(a+b)/2;gc=f(c);if(gc=0) break;else if(gc*gb<0)a=c;ga=gc;elseb=c;gb=gc;cout<<" 方程的根为:X="<<c<&

42、lt;endl; cout<<" *"<<endl;return 0;7.3二分法解非线性方程算法调试结果图示： 用二分法求 在1,2区间上的根图7-2二分法解非线性方程程序算法调试图8.不动点法解非线性方程8.1 不动点法解非线性方程算法概述将改写成对进行迭代，即 其中判断是否成立，成立则返回，不成立就重复以上步骤8.2 不动点法说明函数g(x)的不动点是指一个实数P，满足P=g(P)。从图形角度分析，函数y=g(x)的不动点是y=g(x)和y=x的交点。求不动点的过程可描述如下：第一步：将x0带入迭代公式，作为第一次迭代结果x1第二步：随着不断的

43、迭代，迭代数值会越来越接近不动点的值x0，程序中变量 的类型小数点后的位数是一定的，所以，随着不断的迭代，会出现在误 差线内的两数，那么，此时的可以近似看做方程的根。8.3不动点法解非线性方程程序代码：#include<iostream>using namespace std;#include<cmath>#define MAX 20#define eps 2.2204e-16double g(double x)return 1-x*x/2+x;main()double PMAX=0,err=0.0,relerr=0.0,tol=0.0,p=0.0,p0=0.0;int

44、 k=0,max1=0,i=0;cout<<"不动点法解非线性方程f(x)=1-x2/2+x"<<endl;cout<<"方程在0，1上有解，初始值为p0=0"<<endl;P0=p0=0;max1=100;tol=0.001; for(k=2;k<=max1;k+) Pk-1=g(Pk-2);err=fabs(Pk-1-Pk-2);relerr=err/(fabs(Pk-1+eps);p=Pk-1;if(err<tol)|(relerr<tol)break;if(k=max1) cout

45、<<"迭代次数超过允许的最大迭代次数！"<<endl;cout<<"不动点的近似值为:"<<p<<endl;cout<<"程序迭代次数为"<<k<<endl;cout<<"近似值的误差为:"<<err<<endl;cout<<"求解不动点近似值的序列:"<<endl;for(i=0;i<k;i+) cout<<"

46、"<<Pi;cout<<endl;8.4不动点法解非线性方程程序调试图示用不动点法求解图8-1不动点法解非线性方程程序算法调试图9.牛顿-拉弗森迭代解非线性方程 9.1牛顿-拉弗森迭代解非线性方程算法说明： 存在数pa,b,满足f(p)=0.如果存在一个数对任意近似值，使得由如下迭代定义的序列 收敛到p： （9-1）其中并被称为牛顿-拉弗森迭代函数。由于f(p)=0，显然g(p)=p.这样通过寻找函数g(x)的不动点，可以实现寻找方程f(x)=0的根牛顿-拉弗森迭代。9.2 牛顿-拉弗森迭代解非线性方程流程图NY图9-1 牛顿-拉弗森迭代解非线性方程流程图9.

47、3牛顿-拉弗森迭代解非线性方程程序代码：#include<iostream>#include<cmath>using namespace std;double f(double x)return x*x*x-3*x*x-x+1; double df(double x)return 3*x*x-6*x-1;int main() int k=0,max1=0; double p0,delta,epsilon,p1,err,relerr,y;cout<<"牛顿拉弗森法解非线性方程f(x)=x2-1"<<endl;cout<&l

48、t;"初始值为p0=4.8"<<endl; p0=4.8,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;y=f(p0);if(err<delta)|(relerr<delta)|(fabs(y)<epsilon)break;cout<<"方程的根的近似值为:"<<p0<<endl

49、;cout<<"方程的根的误差估计为:"<<err<<endl;cout<<"迭代次数为:"<<k<<endl;cout<<"方程在p0点的函数值f(p0):"<<y<<endl;return 0;9.4牛顿-拉弗森迭代解非线性方程算法程序调试：用牛顿-拉弗森迭代解非线性方程图9-2牛顿-拉弗森迭代解非线性方程程序算法调试图10.拉格朗日插值10.1 拉格朗日插值算法说明有n+1个点的次数最高为n项的多项式的构造方法,其中拉格朗

50、日（10-1）的系数多项式对每个固定的k，拉格朗日系数多项式具有性质（10-2）10.2 拉格朗日插值算法程序#include<iostream> using namespace std; double func(double X,int k,double x,int n); int main() double Sn=0; int n; cout<<"请输入点的个数n:" cin>>n; double*x=(double*)malloc(n*sizeof(double); double*y=(double*)malloc(n*sizeof

51、(double); double X; int i; for(i=0;i<n;i+) cout<<"请输入x"<<i+1<<",y"<<i+1<<":"<<endl; cin>>xi>>yi; cout<<"请输入x" cin>>X; for(i=0;i<n;i+) Sn=Sn+func(X,i,x,n)*yi; cout<<"通过拉格朗日插值公式所得的结果：&q

52、uot;<<"当x="<<X<<"时，y="<<Sn<<endl; return 0; double func(double X,int k,double x,int n) int i; double Pn=1; double p; for(i=0;i<n;i+) if(i=k) continue; else p=(X-xi)/(xk-xi); Pn=Pn*p; return Pn; 10.3 拉格朗日插值算法程序调试：求解（1，2），（8，9），（0，0），（101）四个点的拉格朗日插值图10-1程序算法调试图设计体会及今后的改进意见