经典四阶龙格库塔法解一阶微分方程组.docx
- 文档编号:27831758
- 上传时间:2023-07-05
- 格式:DOCX
- 页数:47
- 大小:569.69KB
经典四阶龙格库塔法解一阶微分方程组.docx
《经典四阶龙格库塔法解一阶微分方程组.docx》由会员分享,可在线阅读,更多相关《经典四阶龙格库塔法解一阶微分方程组.docx(47页珍藏版)》请在冰豆网上搜索。
经典四阶龙格库塔法解一阶微分方程组
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经典四阶龙格库塔法解一阶微分方程程序代码:
#include
#include
usingnamespacestd;
voidRK4(double(*f)(doublet,doublex,doubley),double(*g)(doublet,doublex,doubley),doubleinitial[3],doubleresu[3],doubleh)
{
doublef1,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;
}
intmain()
{
doublef(doublet,doublex,doubley);
doubleg(doublet,doublex,doubley);
doubleinitial[3],resu[3];
doublea,b,H;
doublet,step;
inti;
cout<<"输入所求微分方程组的初值t0,x0,y0:
";
cin>>initial[0]>>initial[1]>>initial[2];
cout<<"输入所求微分方程组的微分区间[a,b]:
";
cin>>a>>b;
cout<<"输入所求微分方程组所分解子区间的个数step:
";
cin>>step;
cout< : right)< : fixed)< H=(b-a)/step; cout< for(i=0;i {RK4(f,g,initial,resu,H); cout< initial[0]=resu[0];initial[1]=resu[1];initial[2]=resu[2]; } return(0); } doublef(doublet,doublex,doubley) { doubledx; dx=x+2*y; return(dx); } doubleg(doublet,doublex,doubley) { doubledy; dy=3*x+2*y; return(dy); } 1.4经典四阶龙格库塔法解一阶微分方程程序调试结果图示: 应用所编写程序计算所给例题: 其中初值为 求解区间为[0,0.2]。 计算结果为: 图1-2经典四阶龙格库塔法解一阶微分方程算法程序调试图 2.高斯列主元法解线性方程组 2.1高斯列主元法解线性方程组算法分析 使用伪代码编写高斯消元过程: fork=1ton-1do fori=k+1ton l<=a(i,k)/a(k,k) forj=ktondo a(i,j)<=a(i,j)-l*a(k,j) end%endofforj b(i)<=b(i)-l*b(k) end%endoffori end%endoffork 最后得到A,b可以构成上三角线性方程组 接着使用回代法求解上三角线性方程组 因为高斯消元要求a(k,k)≠0(k=1,2,3……n-1)这就需要对高斯消元过程进行完善,即使用高斯列主元法: 其步骤为: 找主元: 计算 并记录其所在行r, 交换第r行与第k行; 以第k行为工具行处理以下各行,使得从第k列的第k+1行到第n行的元素全部为0; 得到增广矩阵的上三角线性方程组; 使用回代法对上三角线性方程组进行求解 2.2高斯列主元法解线性方程组流程图 图2-1高斯列主元法解线性方程组流程图 2.3高斯列主元法解线性方程组程序代码 #include #include #defineN3 usingnamespacestd; voidmain() {inti,j,k,n,p; floatt,s,m,a[N][N],b[N],x[N]; cout<<"请输入方程组的系数"< for(i=0;i {for(j=0;j cin>>a[i][j];} cout<<"请输入方程组右端的常数项: "< for(i=0;i cin>>b[i]; for(j=0;j {p=j; for(i=j+1;i {if(fabs(a[i][j])>fabs(a[p][j])) p=i;} if(p! =j)//交换第p行与第j行 {for(k=0;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 {m=-a[i][j]/a[j][j]; for(n=j;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<<"方程组的解如下: "< for(i=0;i<=N-1;i++) cout< } 2.4高斯列主元法解线性方程组程序调试结果图示: 求解下列方程组 图2-2高斯列主元法解线性方程组程序算法调试图 3.牛顿法解非线性方程组 3.1牛顿法解非线性方程组算法说明 牛顿法主要思想是用 进行迭代。 因此首先需要算出 的雅可比矩阵 ,再求过 求出它的逆 ,当它达到某个精度时即停止迭代。 具体算法如下: 首先设 已知。 : 计算函数 (3-1) : 计算雅可比矩阵 (3-2) (3-3) : 求线性方程组 的解 。 : 计算下一点 重复上述过程。 3.2牛顿法解非线性方程组流程图 图3-1牛顿法解非线性方程组流程图 3.3牛顿法解非线性方程组程序代码 #include #include #defineN2//非线性方程组中方程个数、未知量个数 #defineEpsilon0.0001//差向量1范数的上限 #defineMax100//最大迭代次数 usingnamespacestd; constintN2=2*N; intmain() { voidff(floatxx[N],floatyy[N]);//计算向量函数的因变量向量yy[N] voidffjacobian(floatxx[N],floatyy[N][N]);//计算雅克比矩阵yy[N][N] voidinv_jacobian(floatyy[N][N],floatinv[N][N]);//计算雅克比矩阵的逆矩阵inv voidnewdundiedai(floatx0[N],floatinv[N][N],floaty0[N],floatx1[N]);//由近似解向量x0计算近似解向量x1 floatx0[N]={2.0,0.25},y0[N],jacobian[N][N],invjacobian[N][N],x1[N],errorno rm; inti,j,iter=0; //如果取消对x0的初始化,撤销下面两行的注释符,就可以由键盘向x0读入初始近似解向量 //for(i=0;i //cin>>x0[i]; cout<<"初始近似解向量: "< for(i=0;i cout< cout< do { iter=iter+1; cout<<"第"< //计算向量函数的因变量向量y0 ff(x0,y0); //计算雅克比矩阵jacobian ffjacobian(x0,jacobian); //计算雅克比矩阵的逆矩阵invjacobian inv_jacobian(jacobian,invjacobian); //由近似解向量x0计算近似解向量x1 newdundiedai(x0,invjacobian,y0,x1); //计算差向量的1范数errornorm errornorm=0; for(i=0;i errornorm=errornorm+fabs(x1[i]-x0[i]); if(errornorm for(i=0;i x0[i]=x1[i]; }while(iter return0; } voidff(floatxx[N],floatyy[N]) {floatx,y; inti; x=xx[0]; y=xx[1]; yy[0]=x*x-2*x-y+0.5; yy[1]=x*x+4*y*y-4; cout<<"向量函数的因变量向量是: "< for(i=0;i cout< cout< cout< } voidffjacobian(floatxx[N],floatyy[N][N]) { floatx,y; inti,j; x=xx[0]; y=xx[1]; //jacobianhaven*nelement yy[0][0]=2*x-2; yy[0][1]=-1; yy[1][0]=2*x; yy[1][1]=8*y; cout<<"雅克比矩阵是: "< for(i=0;i {for(j=0;j cout< cout< } cout< } voidinv_jacobian(floatyy[N][N],floatinv[N][N]) {floataug[N][N2],L; inti,j,k; cout<<"开始计算雅克比矩阵的逆矩阵: "< for(i=0;i {for(j=0;j aug[i][j]=yy[i][j]; for(j=N;j if(j==i+N)aug[i][j]=1; elseaug[i][j]=0; } for(i=0;i {for(j=0;j cout< cout< } cout< for(i=0;i { for(k=i+1;k {L=-aug[k][i]/aug[i][i]; for(j=i;j aug[k][j]=aug[k][j]+L*aug[i][j]; } } for(i=0;i {for(j=0;j cout< cout< } cout< 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 {for(j=0;j cout< cout< } cout< 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 {for(j=0;j cout< cout< for(j=N;j inv[i][j-N]=aug[i][j]; } cout< cout<<"雅克比矩阵的逆矩阵: "< for(i=0;i {for(j=0;j cout< cout< } cout< } voidnewdundiedai(floatx0[N],floatinv[N][N],floaty0[N],floatx1[N]) { inti,j; floatsum=0; for(i=0;i {sum=0; for(j=0;j sum=sum+inv[i][j]*y0[j]; x1[i]=x0[i]-sum; } cout<<"近似解向量: "< for(i=0;i cout< cout< 3.4牛顿法解非线性方程组程序调试结果图示: 图3-2牛顿法解非线性方程组程序算法调试图 图3-3牛顿法解非线性方程组程序算法调试图 图3-4牛顿法解非线性方程组程序算法调试图 4.龙贝格求积分算法 4.1龙贝格求积分算法分析 生成j>=k的近似积分结果逼近表,并以R(j+1,j+1)为最终解来逼近积分。 R(j,0)=T(j),j>=0,T(j)为区间逐次减半递推梯形求积分公式算出的结果; R(j,1)=S(j),j>=1,S(j)为区间逐次减半递推辛普森求积分公式算出的结果; R(j,2)=B(j),j>=2,B(j)为递推布尔求积分公式算出的结果; (4-2) 生成 的逼近表 ,并以 为最终解来逼近积分 (4-3) 逼近 存在于一个特别的下三角矩阵中,第0列元素 用基于 个[a,b]子区间的连续梯形方法计算,然后利用龙贝格公式计算 。 当 时,第 行的元素为 当 时,程序在第 行结束。 4.2龙贝格积分算法流程图 图4-1龙贝格积分算法流程图 4.3龙贝格积分算法程序代码 #include usingnamespacestd; #include #include #definef(x)(sin(x))//列举函数 #defineN_H20 #defineMAX10//最大迭代次数 #definea0//所求积分的上下限 #defineb1 #defineepsilon0.0001//所需精度 doubleRomberg(doubleaa,doublebb,longintn) { inti; doublesum,h=(bb-aa)/n;sum=0; for(i=1;i sum+=f(aa+i*h); sum+=(f(aa)+f(bb))/2; return(h*sum); } voidmain() { inti; longintn=N_H,m=0; doubleT[2][MAX+1]; T[1][0]=Romberg(a,b,n); n*=2; for(m=1;m { for(i=0;i {T[0][i]=T[1][i];} T[1][0]=Romberg(a,b,n); n=n*2; for(i=1;i<=m;i++) T[1][i]=T[1][i-1]+(T[1][i-1]-T[0][i-1])/(pow(2,2*m)-1); if((T[0][m-1]-T[1][m]) cout<<"T="< return; } } 4.4龙贝格积分算法调试图 求 在 区间上的精度为0.0001的积分 图4-2龙贝格积分程序算法调试图 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 #defineMAX4 double*diff(doubleX[],intn) {inti=0; double*H=NULL; H=(double*)malloc((n-1)*sizeof(double)); for(i=1;i<=n-1;i++) {H[i-1]=X[i]-X[i-1];} returnH;} double*divide(doubleY[],intN,doubleH[]) {inti=0; double*D=NULL; D=(double*)malloc(N*sizeof(double)); for(i=0;i {D[i]=Y[i]/H[i];} returnD;} main() {double X[MAX]={0,1,2,3},Y[MAX]={0,0.5,2.0,1.5},S[MAX][MAX]={0},temp=0.0,M[MAX]={0};intN=MAX-1,i=0,k=0; doubleA[MAX-1-2]={0},B[MAX-1-1]={0},C[MAX-1-1]={0}; doubledx0=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的三次压紧样条曲线\n\n"); H=diff(X,MAX); D=divide(diff(Y,MAX),N,H); for(i=1;i A[i]=H[i+1]; for(i=0;i B[i]=2*(H[i]+H[i+1]); for(i=1;i C[i]=H[i+1]; U=diff(D,N); for(i=0;i U[i]=U[i]*6; B[0]=B[0]-H[0]/2; U[0]=U[0]-3*(D[0]-dx0); B[N-2]=B[N-2]-H[N-1]/2; U[N-2]=U[N-2]-3*(dxn-D[N-1]); for(k=2;k<=N-1;k++) {temp=A[k-2]/B[k-2]; B[k-1]=B[k-1]-temp*C[k-2]; U[k-1]=U[k-1]-temp*U[k-2];} M[N-1]=U[N-2]/B[N-2]; for(k=N-2;k>=1;k--) M[k]=(U[k-1]-C[k-1]*M[k+1])/B[k-1]
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 经典 四阶龙格库塔法解 一阶 微分 方程组
![提示](https://static.bdocx.com/images/bang_tan.gif)