1.经典四阶龙格库塔法解一阶微分方程组1.1运用四阶龙格库塔法解一阶微分方程组算法分析),,(1k k k y x t f f =， )2,2,2(112g hy f h x h t f f k k k +++=)2,2,2(223g hy f h x h t f f k k k +++=),,(334hg y hf x h t f f k k k +++=),,(1k k k y x t g g =)2,2,2(112g hy f h x h t g g k k k +++=)2,2,2(223g hy f h x h t g g k k k +++=),,(334hg y hf x h t g g k k k +++=)22(6)22(64321143211g g g g hy y f f f f hx x k k k k ++++=++++=++1k k t t h +=+经过循环计算由 推得 ……每个龙格-库塔方法都是由一个合适的泰勒方法推导而来，使得其最终全局误差为()N O h ,一种折中方法是每次进行若干次函数求值，从而省去高阶导数计算。

4阶龙格-库塔方法(RK4)是最常用的，它适用于一般的应用，因为它非常精准，稳定，且易于编程。

000,,t x y ()()111222,,,,t x y t x y （1-1）（1-2）（1-3）（1-4）（1-5）（1-6）（1-7）（1-8）（1-9）（1-10）1.2经典四阶龙格库塔法解一阶微分方程流程图图1-1 经典四阶龙格库塔法解一阶微分方程流程图1.3经典四阶龙格库塔法解一阶微分方程程序代码：#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 initial[3], double resu[3],double h) {double f1,f2,f3,f4,g1,g2,g3,g4,t0,x0,y0,x1,y1;t0=initial[0];x0=initial[1];y0=initial[2];f1=f(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; resu[0]=t0+h;resu[1]=x1;resu[2]=y1;}int main(){double f(double t,double x, double y);double g(double t,double x, double y);double initial[3],resu[3];double a,b,H;double t,step;int i;cout<<"输入所求微分方程组的初值t0,x0,y0:";cin>>initial[0]>>initial[1]>>initial[2];cout<<"输入所求微分方程组的微分区间[a,b]:";cin>>a>>b;cout<<"输入所求微分方程组所分解子区间的个数step:";cin>>step;cout<<setiosflags(ios::right)<<setiosflags(ios::fixed)<<setprecision( 10);H=(b-a)/step;cout<< initial[0]<<setw(18)<<initial[1]<<setw(18)<<initial[2]<<endl; for(i=0;i<step;i++){ RK4( f,g ,initial, resu,H);cout<<resu[0]<<setw(20)<<resu[1]<<setw(20)<<resu[2]<<endl;initial[0]=resu[0];initial[1]=resu[1];initial[2]=resu[2];}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经典四阶龙格库塔法解一阶微分方程程序调试结果图示：应用所编写程序计算所给例题：其中初值为求解区间为[0,0.2]。

计算结果为:图1-2 经典四阶龙格库塔法解一阶微分方程算法程序调试图⎪⎪⎩⎪⎪⎨⎧+=+=yx dtdy y x dt dx232⎩⎨⎧==4)0(6)0(y x2.高斯列主元法解线性方程组2.1高斯列主元法解线性方程组算法分析使用伪代码编写高斯消元过程：for k=1 to n-1 do for 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 j b(i)<=b(i)-l*b(k) end %end of for i end %end of for k最后得到A ，b 可以构成上三角线性方程组 接着使用回代法求解上三角线性方程组因为高斯消元要求a （k,k ）≠0（k=1,2,3……n-1）这就需要对高斯消元过程进行完善，即使用高斯列主元法：其步骤为：①找主元：计算 (){}1max ||()k pk a p k n -≤≤,并记录其所在行r ,||rk a = (){}1max ||()k pk a p k n -≤≤②交换第r 行与第k 行；③以第k 行为工具行处理以下各行，使得从第k 列的第k+1行到第n 行的元素全部为0；④得到增广矩阵的上三角线性方程组； ⑤使用回代法对上三角线性方程组进行求解2.2高斯列主元法解线性方程组流程图图2-1 高斯列主元法解线性方程组流程图2.3高斯列主元法解线性方程组程序代码#include<iostream>#include <cmath>#define N 3using namespace std;void main(){int i,j,k,n,p;float t,s,m,a[N][N],b[N],x[N];cout<<"请输入方程组的系数"<<endl;for(i=0;i<N;i++){for(j=0;j<N;j++)cin>>a[i][j];}cout<<"请输入方程组右端的常数项:"<<endl;for(i=0;i<N;i++)cin>>b[i];for(j=0;j<N-1;j++){ p=j;for(i=j+1;i<N;i++) //寻找系数矩阵每列的最大值{if(fabs(a[i][j])>fabs(a[p][j]))p=i;}if(p!=j) //交换第p行与第j行{for(k=0;k<N;k++){t=a[j][k];a[j][k]=a[p][k];a[p][k]=t;} //交换常数项的第p行与第j行t=b[p];b[p]=b[j];b[j]=t;} //把系数矩阵第j列a[j][j]下面的元素变为0 for(i=j+1;i<N;i++){ m=-a[i][j]/a[j][j];for(n=j;n<N;n++)a[i][n]=a[i][n]+a[j][n]*m; b[i]=b[i]+b[j]*m;}} //求方程组的一个解x[N-1]=b[N-1]/a[N-1][N-1]; //回代法求方程组其他解 for(i=N-2;i>=0;i--) { s=0.0;for(j=N-1;j>i;j--) { s=s+a[i][j]*x[j]; x[i]=(b[i]-s)/a[i][i];}}cout<<"方程组的解如下:"<<endl;for(i=0;i<=N-1;i++) cout<<x[i]<<endl;}2.4高斯列主元法解线性方程组程序调试结果图示：求解下列方程组341-11182522321321321=++=-+=++x x x x x x x x x图2-2 高斯列主元法解线性方程组程序算法调试图3.牛顿法解非线性方程组3.1牛顿法解非线性方程组算法说明牛顿法主要思想是用(1)()()1()()()k k k k x x F x F x +-'=- (0,1,...).k = 进行迭代 。

因此首先需要算出()F x 的雅可比矩阵()F x '，再求过()F x '求出它的逆1()F x -'，当它达到某个精度时即停止迭代。

具体算法如下：首先设k P 已知。

①：计算函数12(,)()(,)k k k k k f p q F P f p q ⎡⎤=⎢⎥⎣⎦②：计算雅可比矩阵1122(,)(,)()(,)(,)k k k k k k k k k f p q f p q x yJ P f p q f p q x y ∂∂⎡⎤⎢⎥∂∂⎢⎥=∂∂⎢⎥⎢⎥∂∂⎣⎦③：求线性方程组 ()()k k J P P F P ∆=-的解P ∆。

④：计算下一点1k k P P P +=+∆重复上述过程。

（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 100 //最大迭代次数using namespace std;const int N2=2*N;int main(){void ff(float xx[N],float yy[N]);//计算向量函数的因变量向量yy[N]void ffjacobian(float xx[N],float yy[N][N]);//计算雅克比矩阵yy[N][N] void inv_jacobian(float yy[N][N],float inv[N][N]);//计算雅克比矩阵的逆矩阵invvoid newdundiedai(float x0[N], float inv[N][N],float y0[N],float x1[N]);//由近似解向量 x0 计算近似解向量 x1floatx0[N]={2.0,0.25},y0[N],jacobian[N][N],invjacobian[N][N],x1[N],errornorm;int i,j,iter=0;//如果取消对x0的初始化，撤销下面两行的注释符，就可以由键盘向x0读入初始近似解向量//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;//计算向量函数的因变量向量 y0ff(x0,y0);//计算雅克比矩阵 jacobianffjacobian(x0,jacobian);//计算雅克比矩阵的逆矩阵 invjacobianinv_jacobian(jacobian,invjacobian);//由近似解向量 x0 计算近似解向量 x1newdundiedai(x0, invjacobian,y0,x1);//计算差向量的1范数errornormerrornorm=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);return 0;}void ff(float xx[N],float yy[N]){float x,y;int i;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;}void ffjacobian(float xx[N],float yy[N][N]){float x,y;int i,j;x=xx[0];y=xx[1];//jacobian have n*n elementyy[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;}void inv_jacobian(float yy[N][N],float inv[N][N]) {float aug[N][N2],L;int i,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;else aug[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;}void newdundiedai(float x0[N], float inv[N][N],float y0[N],float x1[N]) {int i,j;float sum=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]<<" ";cout<<endl;cout<<endl;}3.4牛顿法解非线性方程组程序调试结果图示：图3-2 牛顿法解非线性方程组程序算法调试图图3-3 牛顿法解非线性方程组程序算法调试图图3-4 牛顿法解非线性方程组程序算法调试图4.龙贝格求积分算法4.1龙贝格求积分算法分析生成j>=k 的近似积分结果逼近表，并以R(j+1,j+1)为最终解来逼近积分。