数值分析报告常微分方程的数值解法.docVIP

wilyes11收集 博客(与学习无关)： HYPERLINK /u/1810231802 /u/1810231802 常微分方程的数值解法 实验目的： 目的与要求：通过实验，熟悉常微分方程的数值解法的基本原理。掌握向前欧拉法、向后欧拉法、梯形法、改进欧拉法及三阶、四阶龙格-库塔法等基本算法。 实验内容： 在下列方法中：向前欧拉法、向后欧拉法、梯形法、改进欧拉法及三阶、四阶龙格-库塔法 选择不同的三种算法求下面方程的数值解，并与准确解比较。 以上步长取0.1. 实验要求： 1．编程实现。2.可以自拟实验题，要求同上。 实验程序及结果： 求原函数 DSolve[{y[x]==x^2+x-y[x],y[0]==0},y[x],x] 求得： y=e-x-1+ex-exx+exx2 向前欧拉法 x[0]=0;y[0]=0;h=0.1; f[x_]:=e-x-1+ex-exx+exx2; x[n_]:=n*h; f[u_,v_]:=u2+u-v; K1[n_]:=f[x[n-1],y[n-1]]; y[n_]:=y[n-1]+h*K1[n]; Table[{n,x[n],y[n],f[x[n]],f[x[n]]-y[n]},{n,0,10}]//N; Print[\nn,x[n],y[n],f[x[n]],f[x[n]]-y[n]\n, ,MatrixForm[Table[{n,x[n],y[n],f[x[n]],f[x[n]]-y[n]},{n,0,10}]//N]] Data1=Table[{x[n],f[x[n]]},{n,0,10}]//N; Data2=Table[{x[n],y[n]},{n,0,10}]//N; A0=Plot[f[x],{x,0,1},PlotStyle?{Red,Thick}]; A1=Graphics[{PointSize[0.02],Red,Point[Data1]}]; A2=ListLinePlot[Data2,PlotStyle?{Red,Dashed}]; A3=Graphics[{PointSize[0.02],Blue,Point[Data2]}]; Show[A0,A1,A2,A3,PlotRange?All]  改进欧拉法 f[x_,y_]:=x2+x-y; f[x_]:=e-x-1+ex-exx+exx2; xy:={{0,1}}; h=0.1; Do[xn=xy[[n]][[1]];yn=xy[[n]][[2]]; k1=f[xn,yn];k2=f[xn+h,yn+h*k1]; d=(k1+k2)*h/2 ; Print[n, ,x,n,=,xn+h, ,y,n,=,yn+d, ,f(x,n,)=,f[xn+h], ,f(x,n,)-y(x,n,)=,f[xn+h]-(yn+d)]; xy=Append[xy,{xn+h,yn+d}],{n,1,10}]; A1=Graphics[{PointSize[0.02],Red,Point[xy]}]; A2=ListLinePlot[xy,PlotStyle?Dashed]; A3=Plot[f[x],{x,0,1},PlotStyle?{Red,Thick}]; Show[A3,A2,A1]  四阶龙格-库塔法 f[x_,y_]:=x2+x-y; f[x_]:=e-x-1+ex-exx+exx2;; y0=1;a=0;b=1;n=10; h=(b-a)/n; xx=Table[a+(i-1)h,{i,1,n+1}]//N; y=Table[0,{i,1,n+1}]; y[[1]]=y0//N; For[i=2,i?n+1,i++, a1=h f[xx[[i-1]],y[[i-1]]]; a2=h f[xx[[i-1]]+h/2,y[[i-1]]+a1/2]; a3=h f[xx[[i-1]]+h/2,y[[i-1]]+a2/2]; a4=h f[xx[[i-1]]+h,y[[i-1]]+a3]; y[[i]]=y[[i-1]]+1/6 (a1+2a2+2a3+a4); Print[i-1, ,x,i-1,=,xx[[i]], ,y,i-1,=,y[[i]], ,f(x,i-1,)=,f[xx[[i]]], ,f(x,n,)-y(x,n,)=,f[xx[[i]]]-y[[i]]];] data=Table[{xx[[i]],y[[i]]},{i,1,n+1}]; A1=Graphics[{PointSize[0.02],Blue,Point[data]}]; A2=ListLinePlot[data]; A3=Plot[f[x],{x,0,1},PlotStyle?{Red,Thick}]; Sh