经典四阶龙格库塔法解一阶微分方程组Word文档格式.docx

1、 y1=y0+h*（g1+2*g2+2*g3+g4）/6;resu0=t0+h;resu1=x1;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;coutinitial0initial1initial2;输入所求微分方程组的微分区间a,b:ab;输入所求微分方程组所分解子区间的个数step:step;setiosflags（ios:right）f

2、ixed）setprecision（10）; H=（b-a）/step; initial0setw（18）initial1initial2endl;for（i=0;ii+） RK4（ f,g ,initial, resu,H）; coutresu0setw（20）resu1resu2 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

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

4、end of for i %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高斯列主元法解线性方程组流程图图2-1 高斯列主元法解线性方程组流程图2.3高斯列主元法解线性方程组程序代码#include#define N 3

5、void main（）int i,j,k,n,p;float t,s,m,aNN,bN,xN;请输入方程组的系数N;for（j=0;jbi;for（j=0;N-1; p=j; for（i=j+1;i+） /寻找系数矩阵每列的最大值 if（fabs（aij）fabs（apj） p=i; if（p!=j） /交换第p行与第j行 for（k=0;kk+） t=ajk; ajk=apk;apk=t; /交换常数项的第p行与第j行 t=bp; bp=bj; bj=t; /把系数矩阵第j列ajj下面的元素变为0 m=-aij/ajj; for（n=j;n=0;i-） s=0.0; for（j=N-1;ji

6、;j-） s=s+aij*xj; xi=（bi-s）/aii;方程组的解如下: for（i=0;=N-1;xi2.4高斯列主元法解线性方程组程序调试结果图示：求解下列方程组图2-2 高斯列主元法解线性方程组程序算法调试图 3.牛顿法解非线性方程组3.1牛顿法解非线性方程组算法说明 牛顿法主要思想是用 进行迭代 。因此首先需要算出的雅可比矩阵，再求过求出它的逆，当它达到某个精度时即停止迭代。 具体算法如下：首先设已知。：计算函数（3-1）：计算雅可比矩阵 （3-2）（3-3）：求线性方程组的解。：计算下一点 重复上述过程。3.2牛顿法解非线性方程组流程图图3-1 牛顿法解非线性方程组流程图3.3

7、牛顿法解非线性方程组程序代码#define N 2 / 非线性方程组中方程个数、未知量个数 #define Epsilon 0.0001 / 差向量1范数的上限#define Max 100 /最大迭代次数const int N2=2*N;void ff（float xxN,float yyN）;/计算向量函数的因变量向量yyNvoid ffjacobian（float xxN,float yyNN）;/计算雅克比矩阵yyNNvoid inv_jacobian（float yyNN,float invNN）;/计算雅克比矩阵的逆矩阵invvoid newdundiedai（float x0N,

8、 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读入初始近似解向量/for（ i=0;/x0i;初始近似解向量： for （i=0;x0i do iter=iter+1;第 iter 次迭代开始/计算向量函数的因变量向量 y0ff（x0,y0）;/计算雅克比矩阵 jacobianffjacobian（x0,jacob

9、ian）;/计算雅克比矩阵的逆矩阵 invjacobianinv_jacobian（jacobian,invjacobian）;newdundiedai（x0, invjacobian,y0,x1）;/计算差向量的1范数errornorm errornorm=0; errornorm=errornorm+fabs（x1i-x0i）; if （errornormEpsilon） break; x0i=x1i; while （iterMax）;return 0;void ff（float xxN,float yyN）float x,y; x=xx0; y=xx1; yy0=x*x-2*x-y+0.

10、5; yy1=x*x+4*y*y-4;向量函数的因变量向量是： for（ i=0;yyivoid ffjacobian（float xxN,float yyNN） float x,y; int i,j; /jacobian have n*n element yy00=2*x-2; yy01=-1; yy10=2*x; yy11=8*y;雅克比矩阵是： for（j=0;yyijvoid inv_jacobian（float yyNN,float invNN）float augNN2,L;int i,j,k;开始计算雅克比矩阵的逆矩阵 ：for （i=0; for（j=0; augij=yyij;

11、 for（j=N;N2; if（j=i+N） augij=1; else augij=0;augijk-） for（j=N2-1; augij=augij/augii; invij-N=augij;雅克比矩阵的逆矩阵：invijvoid newdundiedai（float x0N, float invNN,float y0N,float x1N） float sum=0; sum=0; sum=sum+invij*y0j; x1i=x0i-sum;近似解向量：x1i=k 的近似积分结果逼近表，并以R（j+1,j+1）为最终解来逼近积分。R（j,0）=T（j）, j=0,T（j）为区间逐次减半

12、递推梯形求积分公式算出的结果；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龙贝格积分算法程序代码stdio.hmath.h#define f（x） （sin（x） /列举函数#define N_H 20#def

13、ine MAX 10#define a 0 /所求积分的上下限 b 1 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;n; sum+=f（aa+i*h）; sum+=（f（aa）+f（bb）/2; return （h*sum）; long int n=N_H,m=0; double T2MAX+1; T10=Romberg（a,b,n）; n*=2; for（m=1;mMAX;m+）m; T0i=T1i; n=n*2; for （i=1;=m; T1i=T1i-1+（T1i-1-T0i-1）/（pow（2,2*m）-1）; if（T0m-1-T1m）epsilon）T=T1m#define MAX 4 double *diff（double X,int n） int i=0;d