经典四阶龙格库塔法解一阶微分方程组.docx
《经典四阶龙格库塔法解一阶微分方程组.docx》由会员分享，可在线阅读，更多相关《经典四阶龙格库塔法解一阶微分方程组.docx（38页珍藏版）》请在冰点文库上搜索。
经典四阶龙格库塔法解一阶微分方程组
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;jcin>>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范数的上限
#defineMax 100 //最大迭代次数
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;icout< 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（iterreturn0;
}
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;
else aug[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）为区间逐次减半递推辛普森求积分公式算出的结果；
（4-1）
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
#defineMAX 10 //最大迭代次数
#define a 0 //所求积分的上下限
#define b 1
#define epsilon 0.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;
d
升级会员 