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

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

3、=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=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)initial2endl;for

5、(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 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) end %end of for i end %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高斯列主元法解线性方程组流程图图2-1 高斯列主元法解线性方程组流程图2.3高斯列主元法解线性方程组程序代码#include#include #define N 3using namespace std;void main()int i,j

8、,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下面的元素变为0 for(i=j+1;iN;i+) m=-aij/ajj; for

9、(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牛顿法解非线性方程组算法说明 牛顿法主要思想是用 进行迭代 。因此首先需要算出的雅可比矩阵，再求过求出它的逆，当它达到某个精度时即停止迭代。 具体算法如下：首先设已知。：计算函数 （3-1） ：计算雅可比矩阵 （3-2） （3-3）：求

10、线性方程组 的解。：计算下一点 重复上述过程。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);/计算向量函数的因变量向量yyNvoid ffjacobian(float xxN,flo

11、at 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读入初始近似解向量/for( i=0;ix0i;cout初

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

13、rnorm=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 ; coutendl; coutend

14、l;void ffjacobian(float xxN,float yyNN) float x,y; int i,j; x=xx0; y=xx1; /jacobian have n*n element yy00=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 ; cout

16、endl; 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

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

19、ng namespace 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

20、*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;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 MAX 4 double *diff(double X,int n) int i=0; d