四阶龙格库塔法解一阶微分方程组Word文件下载.docx
《四阶龙格库塔法解一阶微分方程组Word文件下载.docx》由会员分享，可在线阅读，更多相关《四阶龙格库塔法解一阶微分方程组Word文件下载.docx（39页珍藏版）》请在冰豆网上搜索。
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];
输入所求微分方程组的微分区间[a,b]:
a>
b;
输入所求微分方程组所分解子区间的个数step:
step;
setiosflags（ios:
:
right）<
fixed）<
setprecision（10）;
H=（b-a）/step;
initial[0]<
setw（18）<
initial[1]<
initial[2]<
endl;
for（i=0;
i<
i++）
{RK4（f,g,initial,resu,H）;
resu[0]<
setw（20）<
resu[1]<
resu[2]<
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<
cmath>
#defineN3
voidmain（）
{inti,j,k,n,p;
floatt,s,m,a[N][N],b[N],x[N];
cout<
请输入方程组的系数"
for（i=0;
N;
{for（j=0;
j<
j++）
cin>
a[i][j];
请输入方程组右端的常数项:
b[i];
for（j=0;
N-1;
{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<
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<
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];
方程组的解如下:
for（i=0;
=N-1;
cout<
x[i]<
2.4高斯列主元法解线性方程组程序调试结果图示：
求解下列方程组
图2-2高斯列主元法解线性方程组程序算法调试图
3.牛顿法解非线性方程组
3.1牛顿法解非线性方程组算法说明
牛顿法主要思想是用进行迭代。
因此首先需要算出的雅可比矩阵，再求过求出它的逆，当它达到某个精度时即停止迭代。
具体算法如下：
首先设已知。
：
计算函数
（3-1）
计算雅可比矩阵
（3-2）
（3-3）
求线性方程组的解。
计算下一点
重复上述过程。
3.2牛顿法解非线性方程组流程图
图3-1牛顿法解非线性方程组流程图
3.3牛顿法解非线性方程组程序代码
#defineN2//非线性方程组中方程个数、未知量个数
#defineEpsilon0.0001//差向量1范数的上限
#defineMax100//最大迭代次数
constintN2=2*N;
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;
//cin>
x0[i];
初始近似解向量：
for（i=0;
x0[i]<
"
do
iter=iter+1;
第"
iter<
次迭代开始"
//计算向量函数的因变量向量y0
ff（x0,y0）;
//计算雅克比矩阵jacobian
ffjacobian（x0,jacobian）;
//计算雅克比矩阵的逆矩阵invjacobian
inv_jacobian（jacobian,invjacobian）;
newdundiedai（x0,invjacobian,y0,x1）;
//计算差向量的1范数errornorm
errornorm=0;
errornorm=errornorm+fabs（x1[i]-x0[i]）;
if（errornorm<
Epsilon）break;
x0[i]=x1[i];
}while（iter<
Max）;
return0;
voidff（floatxx[N],flo
升级会员 