MACM 316 - Assignment 8 Due Date: December 2nd, at 11pm. You must upload both your code (to Assignment 8 scripts/codes) and your report (to Assignment 8 computing report). The assignment is due at 11:00pm. I have set the due time in Canvas to 11:05pm and if Canvas indicates that you submitted late, you will be given 0 on the assignment. Your computing report must be exactly 1 page. There will be a penalty given if your report is longer than one page. Please read the Guidelines for Assignments rst. Keep in mind that Canvas discussions are open forums. Acknowledge any collaborations and assistance from colleagues/TAs/instructor. A. Computing Assignment - Numerical Solution of Kepler's problem One of science's great achievements was the discovery of Kepler's laws for planetary motion; in particular, that the planets follow closed elliptical orbits around the sun. In this assignment you will compare dierent numerical methods for solving Kepler's problem for the motion of a simple solar system consisting of two planets (the two-body problem). For a system of two planets, we may assume one is xed at the origin with the motion of the other planet being in a 2D plane. Let q(t) = q1 (t) q2 (t) , p(t) = p1 (t) p2 (t) , be the position and momentum vectors of the moving planet, respectively. Kepler's laws give the following ordinary dierential equation for q(t) and p(t): q (t) = p(t), p (t) = (q1 (t)2 1 q(t). + q2 (t)2 )3/2 (1) 1. Write a code that implements Euler's method for this problem for time 0 t T = 200 and stepsize h = 0.0005. Use the initial conditions q1 (0) = 1 e, q2 (0) = 0, p1 (0) = 0, p2 (0) = 1+e , 1e e = 0.6. Plot your output in the q1 -q2 plane, i.e. plot the approximate position of the moving planet at time tn for n = 0, 1, . . . , N , where N = T /h . Briey describe the qualitative behaviour of the numerical solution. 2. The ODE (1) has several conserved quantities, including the angular momentum A(t) and Hamiltonian H(t), dened by A(t) = q1 (t)p2 (t) q2 (t)p1 (t), 1 H(t) = (p1 (t)2 + p2 (t)2 ) 2 1 q1 (t)2 + q2 (t)2 . Compute these quantities for your numerical solution. Does your numerical solution also conserve the angular momentum and Hamiltonian? If not, briey comment on their behaviour for large t. 1 3. For systems such as (1), an alternative to the standard Euler's method is the so-called symplectic Euler method h qn+1 . qn+1 = qn + hpn , pn+1 = pn 2 2 (qn+1,1 + qn+1,2 )3/2 Here qn+1,1 and qn+1,2 are the components of the vector qn+1 . Implement this method and compare it with the standard Euler's method. Describe the behaviour of the numerical solution, and also the angular momentum and Hamiltonian. 2 MACM 316 Assignment 7 Kenny Chetal, kca48@sfu.ca 301100888 Aitken'sMethod Without Aitkens Method -0.459360856413897 -0.456828084750435 Q100 Q100 -0.459360905207118 -0.458053579431806 Q200 Q200 -0.459360910486042 -0.458476760127848 Q300 Q300 -0.459360911810678 -0.458691931429551 Q400 Q400 -0.459360912290856 -0.458822404403239 Q500 Q500 -0.459360912505684 -0.458910042739229 Q600 Q600 -0.459360912615926 -0.458973002410678 Q700 Q700 -0.459360912678248 -0.459020440080522 Q800 Q800 -0.459360912716113 -0.459057476931193 Q900 Q900 -0.459360912740438 -0.459087202347749 Q1000 Q1000 -0.459360912806708 -0.459332530323491 Q10000 Q10000 -0.459360912806765 -0.459346611894649 Q20000 Q20000 -0.459360912806770 -0.459351340175269 Q30000 Q30000 -0.459360912806769 -0.459353713945790 Q40000 Q40000 -0.459360912806768 -0.459355142193831 Q50000 Q50000 -0.459360912806767 -0.459356096373886 Q60000 Q60000 -0.459360912806772 -0.459356779084107 Q70000 Q70000 Figure 1: Data without Aitkens method Figure 2: Data using Aitkens method Bonus Question with Aitkens 0.279506125323260 Q100 0.279506100501260 Q200 0.279506097812933 Q300 0.279506097138568 Q400 0.279506096894226 Q500 0.279506096784959 Q600 0.279506096728910 Q700 0.279506096697236 Q800 0.279506096677999 Q900 0.279506096665645 Q1000 0.279506096632022 Q10000 0.279506096631989 Q20000 0.279506096631985 Q30000 0.279506096631985 Q40000 0.279506096631983 Q50000 0.279506096631983 Q60000 0.279506096631984 Q70000 Figure 3: Data for bonus question a) Figure 1 represents an approximation of I using the subintervals as zeros ofthefunctionwithoutAitken's method. The above data for all figures uses the integral command instead of the quad command which is important to note because both may give slightly different answers for each n value. b) The degree of accuracy for Q1000 is 3 digits. I can justify that as I increased the number of n to Q10000 the 4th digit changed. If we want to consider larger values of n, we can see 4 digits of accuracy with this method. c) When using Aitkens method I tried variations of n that you can see above in figure 2. You can see acceleration of convergence using Aitkens method as we compare the results from figure 1 Q70000 to figure 2 Q100. Immediately we can see the increased speedup of convergence in figure 2. Observations of my data would lead to an answer of 13 digits of accuracy. I justify this by looking at my data since Q20000 to Q70000 has 13 digits that remain consistent despite changing n values. a) ForthebonusquestionIchangedmyzero'scalculationsowecalculatethezero'softhebonusfunction.I continued to use the Aitkens method. Iderivedthezero'scalculationfromthefollowingmath. 1. 2. log ) ( 2 log ) ( 2 = 2 + + = 0 (sub in ai for x) 3. log( ) ( 2 log ) ( )^2+ + = 0 2 4. Let bi = -log(ai) 5. -bi * ( )^2 + + 2 = 0 (multiply by -1) 6. bi * ( )^2 2 = 0 b) At this point I analyzed my derived equation and looked what was different between the fzero equation for the earlier parts of the assignment and the needed equation for the bonus. Then I made the necessary changes to get b=fzero(@(x) x*exp(x)^2-i*pi-pi/2,0). Based on the data I collected for the bonus I would say that I have 14 digits accuracy since Q20000 to Q70000 has 14 digits that remain consistent despite changing n values. Furthermore checking with wolfram alpha gave the answer of 0.279506096631990, which adds a bit more justification that another fairly reliable source got a very close answer to mine. 301254393 YouYang Fang Computing Assignment 7 f(x)=(x^-1)*sin((x^-1)*log(x)) n Q(n) Q_hat 100 -0.456828084750435 -0.459360853057809 200 -0.458053579431806 -0.459360904981172 300 -0.458476760127848 -0.459360910440038 400 -0.458691931429551 -0.459360911795865 500 -0.458822404403239 -0.459360912284717 600 -0.458910042739229 -0.459360912502697 700 -0.458973002410678 -0.459360912614303 800 -0.459020440080522 -0.459360912677292 900 -0.459057476931193 -0.459360912715513 999 -0.459634891462477 -0.459360912873714 1000 -0.459087202347749 -0.459360912740043 f(x)=(x^-1)*cos((x^-2)*log(x)) n Q(n) Q_hat 100 0.279490182761134 0.279477453608860 200 0.279510311587285 0.279510496146914 300 0.279503537115601 0.279506672337744 400 0.279506461132379 0.279505158055350 500 0.279505909967502 0.279507415676889 600 0.279505543425124 0.279506256281989 700 0.279506324038014 0.279505542210845 800 0.279506393805645 0.279505924707582 900 0.279506435388056 0.279506063707510 999 0.279505815392537 0.279505815289401 1000 0.279506382141181 0.279505807783769 2000 0.279506032043862 0.279506086160472 3000 0.279506087809366 0.279506087218726 4000 0.279506088315063 0.279506086475413 5000 0.279506091562580 0.279506091541980 6000 0.279505024171042 0.279505027339007 10000 0.279505169671408 0.279505175052342 100000 0.279504788081353 0.279504788037743 From the data I got for the first function, -0.459 is the same starting from n=800. So I should believe that Q(n) has 3 digits of I that is accurately computed. When calculating Q_hat, computer takes a lot shorter amount of time compare to computing Q(n). I try to find out the digit of accuracy of I on Q_hat, n=900 and n=1000 has the same number on the 10th digit, so I compare the successive term of n=1000, which is n=999. Since the 10th digit have different number on n=900 and n=1000, we can conclude that Q_hat has up to 9 digits of accuracy for I starting from n=500. Interestingly, Q_hat has better approximation on the function compare to Q(n). For the second function, I use the exact same code but the different function. However, base on my code, getting an accuracy of 10 digits require more than 100000 iterations which is extremely slow on my computer. MACM316 Assignment 7 Chong Guo / 301295753 3 digits. n (a) Qn 100 200 300 400 500 600 700 800 900 1000 -0.456779250065744 -0.458040888016786 -0.458471019505471 -0.458688667274635 -0.458820299431395 -0.458908572514738 -0.458971917279467 -0.459019606127254 -0.459056815891911 -0.459086665424413 As can be see from Figure 1, when n increases to 1000 we got a value that the rst three digits after the decimal point are relatively stable (hopefully the after increasing of Qn will not change the rst three digits.) To prove this, I used integral function to calculate the integration of this function from a1000 to a20000 (which is a pretty small value) and get the result about 0.00026, this is nearly the truncation error (the integration from 0 to a1000 ) that we have when n = 1000. Note that even we add this value to Q1000 , the rst three digits still remain unchanged, so now we can make sure that we accurately computed 3 digits of I. Figure 1 n (b) Qn 10000 10001 10002 10003 10004 10005 10006 10007 10008 10009 10010 -0.459332524713409 -0.459389298094829 -0.459332530323491 -0.459389292485860 -0.459332535931346 -0.459389286879119 -0.459332541536974 -0.459389281274602 -0.459332547140378 -0.459389275672310 -0.459332552741558 n Qn n Qn 10000 10001 10002 10003 10004 10005 10006 10007 10008 -0.459360912806708 -0.459360912806848 -0.459360912806708 -0.459360912806848 -0.459360912806708 -0.459360912806848 -0.459360912806708 -0.459360912806848 -0.459360912806708 20000 20001 20002 20003 20004 20005 20006 20007 20008 -0.459360912806765 -0.459360912806783 -0.459360912806765 -0.459360912806783 -0.459360912806765 -0.459360912806783 -0.459360912806765 -0.459360912806783 -0.459360912806765 Figure 3 Figure 4 Figure 2 My answer: -0.4593609128067 (13 digits with AbsTol = 1e 15 in integral function) I rst tried with some small n, then keep increasing it until n is about 10000 and got the result like Figure 2 above, we can see that for even number of n, the Qn keeps increasing and the growth rate becomes smaller, for odd number of n, the Qn keeps decreasing and its descent rate becomes smaller too, so we can know that Qn will nally converge to a value between Q10009 (0.459389275672310) and Q10010 (0.459332552741558), but now our accuracy has only 4 digits. After using Aitkens 2 Method we got result like Figure 3 above and it gives us a better accuracy (12 digits). Finally, I tried n to about 20000 and got the answer 0.4593609128067 with 13 correct digits in serveral minutes (Figure 4 ), which is acceptable. (c) (i) We need to change subdivision points an to t our new function. So we need f (x) = 0, which means: cos(x2 log(x)) = 0 = log(x) = k + (k = 0, 1, 2, 3, . . . ) x2 2 = t = k + e2t 2 (Let : x = et , t > 0) (i = 1, 2, 3, . . . ) = te2t i + = 0 2 2 After getting this equation, I changed the root nding part in my code to: = te2t = i + b=fzero(@(x) x*exp(2*x)-i*pi+pi/2,0) a(i)=exp(-b); (ii) With n increase to about 20000, I got the answer with 13 correct digits: 0.2795060966319