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

《经典四阶龙格库塔法解一阶微分方程组(共45页).doc》由会员分享，可在线阅读，更多相关《经典四阶龙格库塔法解一阶微分方程组(共45页).doc（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 #include 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(t0,x0,y0); g1=g(t0,x0,y0); f2=f

3、(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=x1;resu2=y1;int main() double f(

4、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; coutinitial0initial1initial2; coutab; coutstep; coutsetiosflags(ios:right)setiosflags(ios:fixed)setprecision(10); H=(b-a)/step; cout initial0setw(18)initial1setw(18)initial

5、2endl; for(i=0;istep;i+) RK4( f,g ,initial, resu,H); coutresu0setw(20)resu1setw(20)resu2endl; 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 dy;dy=3*x+2*y;return(dy);1.4经典四阶龙格库塔法

6、解一阶微分方程程序调试结果图示： 应用所编写程序计算所给例题： 其中初值为求解区间为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*b(k) end %end of for iend %end of for k最后得到A，b可以

7、构成上三角线性方程组接着使用回代法求解上三角线性方程组 因为高斯消元要求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 高斯列主元法解线性方程组流程图2.3高斯列主元法解线性方程组程序代码#include#include #define N 3u

8、sing namespace std;void main()int i,j,k,n,p;float t,s,m,aNN,bN,xN;cout请输入方程组的系数endl;for(i=0;iN;i+)for(j=0;jaij;cout请输入方程组右端的常数项:endl;for(i=0;ibi;for(j=0;jN-1;j+) p=j; for(i=j+1;ifabs(apj) p=i; if(p!=j) /交换第p行与第j行for(k=0;kN;k+)t=ajk;ajk=apk;apk=t; /交换常数项的第p行与第j行 t=bp; bp=bj; bj=t; /把系数矩阵第j列ajj下面的元素变为

9、0for(i=j+1;iN;i+) m=-aij/ajj; for(n=j;n=0;i-) s=0.0;for(j=N-1;ji;j-)s=s+aij*xj; xi=(bi-s)/aii; cout方程组的解如下:endl;for(i=0;i=N-1;i+) coutxiendl;2.4高斯列主元法解线性方程组程序调试结果图示：求解下列方程组 图2-2 高斯列主元法解线性方程组程序算法调试图 3.牛顿法解非线性方程组3.1牛顿法解非线性方程组算法说明 牛顿法主要思想是用 进行迭代 。因此首先需要算出的雅可比矩阵，再求过求出它的逆，当它达到某个精度时即停止迭代。 具体算法如下：首先设已知。：计算

10、函数 （3-1） ：计算雅可比矩阵 （3-2） （3-3）：求线性方程组的解。：计算下一点重复上述过程。3.2牛顿法解非线性方程组流程图图3-1 牛顿法解非线性方程组流程图3.3牛顿法解非线性方程组程序代码#include#include#define N 2 / 非线性方程组中方程个数、未知量个数 #define Epsilon 0.0001 / 差向量1范数的上限#define Max 100 /最大迭代次数using namespace std;const int N2=2*N;int main()void ff(float xxN,float yyN);/计算向量函数的因变量向量yyN

11、void ffjacobian(float xxN,float yyNN);/计算雅克比矩阵yyNNvoid inv_jacobian(float yyNN,float invNN);/计算雅克比矩阵的逆矩阵invvoid newdundiedai(float x0N, float invNN,float y0N,float x1N);/由近似解向量 x0 计算近似解向量 x1float x0N=2.0,0.25,y0N,jacobianNN,invjacobianNN,x1N,errornorm;int i,j,iter=0;/如果取消对x0的初始化，撤销下面两行的注释符，就可以由键盘向x0

12、读入初始近似解向量/for( i=0;ix0i;cout初始近似解向量：endl;for (i=0;iN;i+)coutx0i ; coutendl;coutendl;doiter=iter+1;cout第 iter 次迭代开始endl;/计算向量函数的因变量向量 y0ff(x0,y0);/计算雅克比矩阵 jacobianffjacobian(x0,jacobian);/计算雅克比矩阵的逆矩阵 invjacobianinv_jacobian(jacobian,invjacobian);/由近似解向量 x0 计算近似解向量 x1newdundiedai(x0, invjacobian,y0,x1

13、);/计算差向量的1范数errornormerrornorm=0;for (i=0;iN;i+)errornorm=errornorm+fabs(x1i-x0i);if (errornormEpsilon) break;for (i=0;iN;i+)x0i=x1i; while (iterMax);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;iN;i+) coutyyi ;

14、coutendl;coutendl;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;iN;i+) for(j=0;jN;j+) coutyyij ;coutendl; coutendl;void inv_jacobian(float yyNN,float invNN)float augNN2,L; int i,j,k; cout

15、开始计算雅克比矩阵的逆矩阵 ：endl; for (i=0;iN;i+) for(j=0;jN;j+) augij=yyij; for(j=N;jN2;j+)if(j=i+N) augij=1;else augij=0; for (i=0;iN;i+) for(j=0;jN2;j+) coutaugij ; coutendl;coutendl;for (i=0;iN;i+) for (k=i+1;kN;k+) L=-augki/augii;for(j=i;jN2;j+) augkj=augkj+L*augij; for (i=0;iN;i+) for(j=0;jN2;j+) coutaugij

16、 ; coutendl;cout0;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;iN;i+) for(j=0;jN2;j+) coutaugij ; coutendl;cout=0;i-) for(j=N2-1;j=0;j-)augij=augij/augii;for (i=0;iN;i+) for(j=0;jN2;j+) coutaugij ; coutendl; for(j=N;jN2;j+) invij-N=augij;coutendl;cout雅克比矩阵的逆

17、矩阵： endl;for (i=0;iN;i+) for(j=0;jN;j+) coutinvij ; coutendl;coutendl;void newdundiedai(float x0N, float invNN,float y0N,float x1N)int i,j;float sum=0;for(i=0;iN;i+) sum=0; for(j=0;jN;j+) sum=sum+invij*y0j;x1i=x0i-sum; cout近似解向量：endl;for (i=0;iN;i+) coutx1i ; coutendl;cout=k 的近似积分结果逼近表，并以R(j+1,j+1)为

18、最终解来逼近积分。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）生成的逼近表，并以为最终解来逼近积分（4-3）逼近存在于一个特别的下三角矩阵中，第0列元素用基于个a,b子区间的连续梯形方法计算，然后利用龙贝格公式计算。当时，第行的元素为当时，程序在第行结束。4.2龙贝格积分算法流程图 图4-1 龙贝格积分算法流程图4.3龙贝格积分算法程序代码#includeusing na

19、mespace std;#include#include#define f(x) (sin(x) /列举函数#define N_H 20#define MAX 10 /最大迭代次数#define a 0 /所求积分的上下限#define 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;in;i+) sum+=f(aa+i*h); sum+=(f(aa)+f(bb)/2; return (h*sum)

20、;void main() int i; long int n=N_H,m=0; double T2MAX+1; T10=Romberg(a,b,n); n*=2; for(m=1;mMAX;m+) for(i=0;im;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) coutT=T1mendl;return; 4.4龙贝格积分算法调试图求在区间上的精度为0.0001的积分 图4-2 龙贝格积分程序算法调试图

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

22、AX 4 double *diff(double X,int 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;iN;i+) Di=Yi/Hi; return D; main() doubleXMAX=0,1,2,

23、3,YMAX=0,0.5,2.0,1.5,SMAXMAX=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

24、N-2;i+) Ai=Hi+1; for(i=0;iN-1;i+) Bi=2*(Hi+Hi+1); for(i=1;iN-1;i+) Ci=Hi+1; U=diff(D,N); for(i=0;iN;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=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

25、=(Mk+1-Mk)/(6*Hk); Sk1=Mk/2; Sk2=Dk-Hk*(2*Mk+Mk+1)/6; Sk3=Yk; printf(求得的三次紧压样条曲线的矩阵S为:n); for(i=0;iMAX-1;i+) for(k=0;kMAX;k+) printf(%lft,Sik); printf(n); for(i=0;i x1=0:.01:1; y1=polyval(S(1,:),x1-X(1);x2=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,y

26、3,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)绝对值的最大值，即误差向量的范数；二是误差绝对值的和，即误差向量r的1范数；三是误差平方和的算术平方根，即误差向量r的2范数；前两种方法简单、自然，但不便于微分运算 ，后一种方法相当于考虑 2范数的平方，因此在曲线拟合中常采用误差平方和来 度量误差

27、(i=0，1，m)的整体大小。数据拟合的具体作法是：对给定数据 (i=0,1,，m)，在取定的函数类中,求,使误差(i=0,1,m)的平方和最小。6.2 M次多项式曲线拟合流程图FT开始输入x0、输出极小值点和极小值结束图6-1 M次曲线多项式曲线拟合流程图6.3 M次多项式曲线拟合算法程序代码#include #include #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;iMAX;i+) for(j=0;jMAX;j+) if(i=

28、j) Eij=1; for(i=0;in-1;i+) temp=Xii; for(j=0;jn;j+) Xij=Xij/temp; Eij=Eij/temp; for(k=0;kn;k+) if(k=i) continue; temp=-Xii*Xki; for(j=0;jn;j+) 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; prin

29、tf(tttM次多项式曲线拟合nn输入拟合的点的个数:); scanf(%d,&n); printf(n输入%d个点的X坐标序列:n,n); for(i=0;in;i+) scanf(%lf,&Xi); printf(n输入%d个点的Y坐标序列:n,n); for(i=0;in;i+) scanf(%lf,&Yi); printf(n输入拟合次数:); scanf(%d,&M); for(i=0;in;i+) for(k=1;k=M+1;k+) Fik-1=pow(Xi,k-1); for(i=0;in;i+) for(j=0;jM+1;j+) BFji=Fij; for(i=0;iM+1;i

30、+) for(j=0;jM+1;j+) for(k=0;kn;k+) Aij+=BFik*Fkj; for(i=0;iM+1;i+) for(j=0;jn;j+) Bi+=BFij*Yj; inv(A,n,E); for(i=0;iM+1;i+) for(j=0;j=0;i-) printf(%lft,Ci); printf(n拟合后的%d次多项式为:n,M); printf(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次多项式曲线拟合算法程

31、序调试结果：求解（-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.二分法解非线性方程7.1 二分法解非线性方程算法说明：.是起始区间，是中点。.是第二个区间，它包含零点r，同时是中点，区间的宽度范围是的一半。.得到第n个区间（包含r，并有中点）后，可构造出，它也包括r，宽度范围是的一半。 7.2二分法解非线性方程程序代码#includeusing namespace std;#define eps

32、 0.001double f(double x) return 2-x*x;int main()double ga,gb,gc,a,b,c; cout *endl;cout 用二分法寻找方程2-x*x=0的根endl;cout 求根区间为(1,2)endl; cout *eps)c=(a+b)/2;gc=f(c);if(gc=0) break;else if(gc*gb0)a=c;ga=gc;elseb=c;gb=gc;cout 方程的根为:X=cendl; 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第二步：随着不断的迭代，迭代数值会越来越接近不动点的值x0，程序中变量