数值计算课程设计
1 经典四阶龙格库塔法解一阶微分方程组
1.1算法说明：
4阶龙格-库塔方法(RK4)可模拟N=4的泰勒方法的精度。这种算法可以描述为，自初始点(t0,y0)开始，利用
yk?1?yk?生成的近似序列，其中
h(f1?2f2?2f3?f4)
6f1?f(tk,yk)
hhf2?f(tk?,yk?f1)22 hhf3?f(tk?,yk?f2)22f4?f(tk?h,yk?hf3)1.2 用龙格库塔法求解求解微分方程
?x?x?2y??y?3x?2y?x?0??6 满足初值条件：?
??y0?4?1.3 算法流程图：
图1-1 四阶龙格库塔法解一阶微分方程组
1
经典四阶龙格库塔法解一阶微分方程组
1.4 程序调试：
组建调试，确保程序可运行时输入初值，区间，步长。 1.5 程序运行运行界面及运行结果：
图1-2 四阶龙格库塔法解一阶微分方程组 1.6 源程序代码: #include
float ff(float t,float xx,float yy) {
xx=xx+2*yy; return xx; }
float gg(float t,float xx,float yy) {
yy=3*xx+2*yy; return yy; }
void rks4(float xx[N],float yy[N],float tt[N],double a,double b) { float h,p1,p2,p3,p4,q1,q2,q3,q4; xx[0]=6; yy[0]=4; int i,p; h=(b-a)/N;
for(p=0;p 2 数值计算课程设计 tt[p]=a+h*p; for(i=0;i { p1=h*ff(tt[i],xx[i],yy[i]); q1=h*gg(tt[i],xx[i],yy[i]); p2=h*ff(tt[i]+h/2,xx[i]+p1/2,yy[i]+q1/2); q2=h*gg(tt[i]+h/2,xx[i]+p1/2,yy[i]+q1/2); p3=h*ff(tt[i]+h/2,xx[i]+p2/2,yy[i]+q2/2); q3=h*gg(tt[i]+h/2,xx[i]+p2/2,yy[i]+q2/2); p4=h*ff(tt[i]+h,xx[i]+p3,yy[i]+q3); q4=h*gg(tt[i]+h,xx[i]+p3,yy[i]+q3); xx[i+1]=xx[i]+(p1+2*p2+2*p3+p4)/6; yy[i+1]=yy[i]+(q1+2*q2+2*q3+q4)/6; } } int main() { float xx[N],yy[N],tt[N]; rks4(xx,yy,tt,0,0.2); int i,j,k; for(i=0;i cout< for(j=0;j cout< for(k=0;k cout< 3 高斯列主元法解线性方程组 2 高斯列主元法解线性方程组 2.1 算法说明： Gauss列主元序消元法主要根据线性方程组人以交换两个方程的次序，方程组的同解解性不变，且解的分量次序也不变。于是，第k步在顺序学消元法进行之前，从Ak的第k’列元素a[k][k],a[k+1][k],……a[n][k]中选出绝对值最大者，并进行记录所在行，即 |a[i][k]|=max|a[i][k]| 如果l不等k，则交换矩阵的第k’行与第l列所对应的元素，然后，再进行第k步顺序消元法。 2.2 算法流程图： 4 数值计算课程设计 图2-1 高斯列主元法解线性方程组 2.3高斯列主元程序调试： 对所编写的高斯列主元程序进行编译和链接，然后执行得如下所示的窗口，我们按命令输入增广矩阵的行数为3，输入3行4列的增广矩阵运行。 2.4 程序运行运行界面及运行结果： 图2-2 高斯列主元法解线性方程组 5