计算物理第六次作业.docVIP

计算物理第六次作业 ——编写程序求解微分方程 问题描述： 已知的一个微分方程及初值条件如下所示： 编写一个程序求解该微分方程，并分别以步长0.1、0.1/2、0.1/4、0.1/8计算y(1.5)的值，然后以图像的形式将通过程序运算得到的结果与精确解进行比较。 主要算法与公式： 在这里，我们求解微分方程的方法是四阶Runge-Kutta法，其公式如下所示： 采用这种方法计算的结果一般比较精确，并且编写的程序具有占用内存小的优点。 源代码： （Runge-Kutta） #includeiostream #includecmath #includeiomanip #include fstream using namespace std; void main() { //ofstream cout(GM1.txt,ios::out); coutsetprecision(8)setiosflags(ios::fixed)setiosflags(ios::right); double S[4]={0,0,0,0}; double n=15; for(int i=0;i4;i++) { double K1,K2,K3,K4; double M,N=3,h=1.5/n; coutsetw(8)xsetw(12)ysetw(12); coutxsetw(12)ysetw(12)xsetw(12)yendl; for(int j=0;jn;j++) { double k=j; K1=-pow((k*h),2)*pow(N,2); K2=-pow((k*h+h/2),2)*pow((N+h*K1/2),2); K3=-pow((k*h+h/2),2)*pow((N+h*K2/2),2); K4=-pow((k*h+h),2)*pow((N+h*K3),2); M=N+h*(K1+2*K2+2*K3+K4)/6; coutsetw(12)k*hsetw(12)N; if(j%3==2) coutendl; N=M; } coutsetw(12)1etw(12)Nendl; S[i]=N; n=n*2; } for(int i=0;i4;i++) { coutsetw(12)S[i]endl; } } 流程图： 运行的结果： 当h=0.1时，运行的结果如下图所示： 观察运行的结果图，当x=1.5时，y的值是0 当h=0.1/2时，运行的结果如图所示： 所以，通过运行的结果，y(1.5)的值为0 当h=0.1/4时，运行的结果如图所示： 从上图中得到y(1.5)的值是0 当h=0.1/8时，运行的结果如图所示： 运行的结果告诉我们y(1.5)=0 综上所述，在不同的步长下，算得y(1.5)结果： h 0.1 0.05 0.025 0.0125 y(1.5) 0000 误差分析与绘图： 已知该微分方程的精确解是，所以，我们将在不同步长通过Runge-Kutta得到的结果与精确值对比，并研究与精确解的误差状况。 (1)当h=0.1时，运行的结果与精确解的关系： 各点的结果与精确解的偏离状况： (2)当h=0.05时，运行的结果与精确解的关系： 各点的结果与精确解的偏离状况： (3)当h=0.025时，运行的结果与精确解的关系： 各点的结果与精确解的偏离状况： (4)当h=0.0125时，运行的结果与精确解的关系： 各点的结果与精确解的偏离状况： 分析：比较了不同步长下的结果，发现结果与精确值得相对偏离都是比较小的，几乎可以忽略，而且取得步长越小，偏离越小，最后只在0附近，得到的结果也越接近精确解。