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
copyright@ 2008-2022 冰豆网网站版权所有
经营许可证编号:鄂ICP备2022015515号-1