數(shù)值大作業(yè)第三次

x(t)y(t)x(0)x'(t)ay(t

x'(t)y'(t)x''(t)abx(t)y''(t)aby(t)r2abrxc1cos(t)d1yc2cos(t)d2x(0)0,y(0)bc10,c2bxd1ybcos(t)d2x'(t)d1cos(t)abcos(t)ad2d1則

ab,d2xybx2y2 利用①改進的歐拉折線法，②4h，的時間步要足夠多――至少10000步）x'f1(x,y)y'f(x,y) xpxnhf1(xn,yn)xnxcxnhf1(xp,yn1)xnxx

)/2

ha(

)/ypynhf2(xn,yn)ynycynhf2(xn1,yp)ynyy

(y

)/2

)/ ((4abh2)x4ahy)/(4abh2

((4abh2)yn4bhxn)/(4abh2xn1xnh(K12K22K3K4)/

ynh(L12L22L3L4)/K1f1(xn,yn)aynLf(x,y) K2f1(xnhK1/2,ynhL1/2)a(ynhL1/ 2

hK/2,yhL/2)

hK1/L K3f1(xnhK2/2,ynhL2/2)a(ynhL2/Lf(xhK/2,yhL/2)b(xhK/

K4f1(xnhK3,ynhL3)a(ynhL3Lf(xhK,yhL3)b(xhK 并設定步長h和步數(shù)，便可根據(jù)上述迭代公式進行迭代和計算。x2y2 %a,b%h%step時間步總數(shù)fory(n+1)=((4-(h^2)*a*b)*y(n)-4*h*b*x(n))/(4+(h^2)*a*b);y1=zeros(step+1,1);forK1=a*y1(n);L1=-b*x1(n);holdon;symsu set(k,'linestyle',':','color','k','linewidth',2);plot(x,y,'r-','linewidth',2);%繪制改進的歐拉折線法計算軌道axis([-55-33]);legend('AccurateEllipse','ImprovedEuler','RungeKutta4'); holdoff;321y0y xyy3.9833.98353.9843.98453.9853.98553.9863.98653.987x4階標準龍格-庫塔法都能保持橢圓軌道不變；321y0y xyxy4階標準龍格-庫塔法都能保持橢圓軌道不變；321y0y xyy x4階標準龍格-庫塔法都能保持橢圓軌道不變；321y0y xyy x4階標準龍格-庫塔法都能保持橢圓軌道不變；321y0y xyyyy

3.8623.8643.8663.8683.873.8723.8743.8763.878x x4階標準龍格-庫塔法隨著迭代的進行越來越偏4階標準龍格-庫塔法計算的軌道321y0y x1yy0 x321y0y x4階標準龍格-庫塔法計算時，原方步長較小時，隨著迭代的進行，改進的歐拉折線法和4階標準龍格-庫塔法①步長相同時，不同方法對橢圓軌道的影響：隨著迭代進行，改進的歐拉折②方法相同時，不同步長對