華南理工大學《計算方法》實驗報告

1、計算機科學與工程學院專業(yè)計算機科學與技術(shù)（全英創(chuàng)新班）課程編號145022#諭理乂人第south china university of technology華南理工大學計算方法實驗報告 華南理工大學計算方法實驗報告課程學分3分起始日期2015年5月28日實驗題目一solution of nonlinear equationsbracketing method and newton-raphson method姓名:xxxxxx學號：xxxxxxxxxxxx時間：2015年5月28日【題目】1. find an approximation (accurate to 10 decimal pla

2、ces) for the interest rate i that will yield a total annuity value of $500,000 if 240 monthly payments of $300 are made.2. consider a spherical ball of radius r = 15 cm that is constructed from a variety of white oak that has a density of p = 0.710. how much of the ball (accurate to eight decimal pl

3、aces) will be submerged when it is placed in water?3. to approximate the fixed points (if any) of each function. answers should be accurate to 12 decimal places produce a graph of each function and the line y = x that clearly shows any fixed points.(a) g(x) = x5 - 3x3 一 2x2 + 2(b) g(x) = cos(sin(x)【

4、實驗分析】1. assume that the value of the k-th month is ak. the value next month isak+1 = 300 + a* * 人 where a0 = 300so the total annuity value after 240 months is :(1 -(1 + /)240)k=239a =ak = 300 + 300 * (1 + /) + + 300 * (1 + i)239 = 300 *k=0to find the rate i that satisfied a=5000000, we can find the

5、solution to theequation g(i) = a- 500000 = 300 *(打：穿)一 500000 = 0the bisection method can be used to find the solution of g(i) = 0.rdmw = i 7r(r2 (% r)2)dx 丿o2. the massmw of water displaced when a sphere is submerged to a depth d is: 7rd2(3r d)3the mass mb of the ball ismb = happlying archimedes5 l

6、aw, mw = mb, producethe following equations:7r(d3 3d2r + 4r3p)3so g(d)=d3 一 3d2r + 4r3p, and bisection method can be used to find the root.3. in this problem, the target is to find the fixed point by fixed point iteration. so we can use fixed point iteration method to solve the problem.(a) g(x) = x5

7、 3x3 2x2 + 2,so g'(x) = 5x4 9x2 4%，g"(x) = 20x3 18% 4plot the function g(x) and y=x, it can be seen that there are three fixed point,where 2 << 1.6,0.6 < x2 < 0.7, 1.9 < x± < 2.1it can be solved that:for2 < x<1.6,g(x) <> 9(一1.6) = 16.128for0.6 <x<0.7

8、,g(x)<0,9©) < /(07) = -6.0095forl.9 <%<2.1,心)>0,9©) > (1.9) = 25.0705soj/crj > 1,0(兀2)1 > 1 ,0(兀 3)i > 1，they are all repelling fixed point, fixed point iteration method cannot converge to them.(b) g(x) = cos (sin (x) from the graph below we know that the fixed p

9、oint lies between 0.5 and 1. where g'(%) = cos(x) sin(sin(x) < 1. so the fixed point is attractive fixed point. i initially guess that p=0.75, and apply fixed point iteration method to find the solution.【實驗過程和程序】1. the value of g(i) will change significantly but little change of i, so that it

10、's unwise to plot the figure about i and g(i). instead, the interest will be bound by 0,1 by common sense. so i set the left end point of bisection method to 0, right end point to 1. the precision is set to le-10.code for g(i) is :function result=g(i) result=500000-300»(1-(1+1/12).*240)/(1-

11、(1+i);code for applying bisection method is :czeropoint, yc, err2"bisect (f 1_1, 0,1, 0. 001): fprintf( the ratio is b 10fn , zeropoint);2. the since the density of the ball is smaller than water, so depth submerged should not be larger than 2r, and not less than 0. so i set the left end point

12、of bisection method to 0 and right to 2r.code for g(d) is :funct ion result=g(d)rou» 0 71;r«0.15; result=4*atan(l. 0)»(d. 3-3*d. 2»r+4*r 3*rou)/3;code for applying bisection method is :zzeropoint.yc, err>bisect <©f 1_2, 0, 2*r, 0. 001);3. (a)code fbr function g:1functi

13、on y=x *5-3*x. '3-2*x. 2+2:code for applying fixed point iteration: 'k, fixpoint, err, pl«f lxpt (©f l_3_a, 2, le-12,100):code for plot the graph:clear, clc, closeh畫盹壯曲紅 xw-2. 5: 0. 01: 2. 6; 嚴;plot (x,y> b', cmin(x>» naz(x)l 2nin(x),bax(x)>' f )： xlabelc x1 )

14、; ylabel (' y*);xlm(max 氏)：) ylin( lmm(y) nax (y) 1);grid on： hold on： %龍,法求紛乘點ik.fixpoint.err.pl=flxpt(©f1.3.a. 2. le-12. 100): plot (cf ixpoint. fixpoint： zmn(y) f ixpointl r ) fprintf (* fixed point is 12fn*. fixpoint):(b)code of this problem is similar to the previous one, so ifs not il

15、lustrated here.【實驗結(jié)果】1. result for the program is :the ratio is 0.3178710938so the annual interest rate is 31.78710938%2. x-axis is the depth and y-axis is for g(d). blue line is the function value versus depth. red line is y=0. it can be seen from the graph that the solution is about 19. and by usi

16、ng bisection method, the solution depth= 19.306641 cm is found, with 8 significance bit.the ball will be subserged for 19.306641cm3. (a)no fixed point can be found. by solving this equation, we can find that x=2 is one of the solution. if we set the initial point to 2.00000000001 or 1.9999999, the r

17、esult is the same:ed point is nan f we set initial point to 2, the point canbe found: fixed point is 2. oouooooooooo this is consistent with the conclusion discussed in the previous part(b)the fixed point p=0.768169156736 is found:fixed point is 0.768169156736實驗題目二solution of linear systemsgaussian

18、elimination and pivoting姓名:xxxxxx學號：xxxxxxxxxxxx時間：2015年5月28 口【題目】find the sixth-degree polynomial y = al + a2x + a3xa2 + a4xa3 + a5xa4 + a6xa5 + a7xa6 that passes through (0, 1), (1, 3), (2, 2), (3, 1), (4, 3), (5, 2), and (6, 1). use the plot command to plot the polynomial and the given points on

19、the same graph. explain any discrepancies in your graph【實驗分析】to find the coefficient of each term, we can substitute all the given point to the equations and construct an argument matrix ax=b. then perform triangular factorization and solve the upper-triangular and lower triangular matrix by back su

20、bstitution method.since det(a)= 24883e+07#), there exist a unique solution for ax=b.【實驗過程和程序】code for construct coefficient matrix a:a=zeros (7, 7):for i=l:7for jl:7a(i, j)-(i-l) (j-1):endendcode for constructmatrix b:2,1,3, 2,1'solve x by triangular factorization and back substitution method:x=

21、uptrbk(a, b):【實驗結(jié)果】the result for coefficient matrix:ans =1. 0000 1.8000 11.0250 -10.56253.9375 -0.63750.0375the graph of the polynomial is shown below. blue lines denotes the graph for thepolynomial and red points are the sample points. this polynomial fits the points well.實驗題目三iteration of linear

22、system solutionseidel iteration姓名:xxxxxx學號：xxxxxxxxxxxx時間：2015年5月28 口【題目】(a) start with p0 = 0 and use jacobi iteration to find pk for k= 1, 2, 3. will jacobi iteration converge to the solution?(b) start with p0 =0 and use gauss-seidel iteration to flndp£ for k = 1, 2, 3. willgauss-seidel itera

23、tion converge to the solution?【實驗分析】in this problem we only need to construct the coefficient matrix and matrix b and then apply the jacobi iteration method and gauss-seidel iteration method to find the solution. since |an|> a12 + |a13|,|a22l> 1211 + m23i，1也31> lsil +a satisfy strictly diag

24、onally dominant, so the iteration steps will converge.【實驗過程和程序】a=5,-l, 1;2,1,4j:b=c10,:0t ;for i-l:3p_k (:, i, l)-jacobi (a, b, p.init, 0. 01, i):p_k (:, i, 2)=gseid(a,b, p.init, 0. 01, i):end【實驗結(jié)果】points in first three iteration is shown below. both this two method will converge.pl for jacobipl for

25、 gauss-seideliterationiteration2,0000.1.3750.0. 75002.0000. 0.8750,1.0313p2 for jacobip2 for gauss-seideliteration iteration2.1250, 0.9688,0. 90631.9688, 1.0117,0. 9893p3 for jacobip3 for gauss-seideliteration iterationis2,0125. 0. 9570»1<03912.0045. 0.9975.1.0017實驗題目四interpolationspline int

26、erpolation姓名:xxxxxx學號：xxxxxxxxxxxx時間：2015年5月28 口【題目】the measured temperatures during a 5-hour period in a suburb of los angeles onnovember 8 are given in the following tabic.(a) to constract a lagrange interpolatory polynomial for the data in the table.time, p.m.degrees fahrenheit166266365464s63663(

27、b) to estimate the average temperature during the given 5-hour period.(c) graph the data in the table and the polynomial from part (a) on the same coordinate system. discuss the possible error that can result from using the polynomial in part (a) to estimate the average temperature.【實驗分析】(a) x=l，2,3

28、,4,5,6 y=66,66,65,64,63,63，the target is to construct a lagrange interpolatory polynomial for x and y. we can apply the lagrange approximation method to solve it.(b) after we get the interpolatory polynomial, we can use a great amount of sample points and estimate the function value of them and take

29、 the mean.(c) error contain in lagrange interpolatory polynomial of 5-degree is:6!e5 (%) =(宀。)(宀】)(宀2)(3)(沙宀5)嚴(0, where xo<c<%5error on the average temperature equals to the mean of e5(%), which is：error = -es (x)dxx5-x0jx071f%5 o 一 兀o)(x 一 衍)(兀一兀2)(兀一兀3)0 一 咒4)(咒 一 x5)f6(c)=站卅"改6!xq嚴(c)

30、(%5 兀 0)6!/6(c) fx6 + 21x5 + 175%4 一 735x3 + 1624* 一 1764% + 720)必 嚴(c) r 7 + zx6 + 35x5 + 空尤4 + 竺織3 + 882* + 720x )r243711=240.06嚴(c)rx5i (% xo)(x %!)(% x2)(x _ 兀3)(無 _ %4)(% _ x5)dx xo5*6!5*6!【實驗過程和程序】(a) construct x and y, and solve it by lagrange approximation:x=l:6;y二66, 66, 65,64,63,63: plotho

31、ld onzc,lz=lagran(x, y):(b) after we get the polynomial function, we calculate function value on a set of sample points and get their mean as the average temperature.x_temp= (1: 0. 01: 6)':x_new=ones(length(x_temp), 1);ifor i»l:5x.,new-x.teap. i.x.new;endy_new= (x_ne"c )':avg_tempe

32、r=bean(y_new)(c) plot x_new and y_new in part(b) is equivalent to plot pn(x). and error bound can only be obtain by analyzing since f(x) is unknown.【實驗結(jié)果】(a) we can get the result:c =0.0167-0.29172.0000-6.70839.983361.0000which means that:f(x) = 0.0167x5 一 0.2917%4 + 2x3 一 6.7083%2 + 9.9833% + 61(b)

33、 by taking the mean, the average temperatureis 64.5°f:avg_te«per =64 5000(c) the graph for polynomial and the sample points is shown below. error bound forthe average temperature is 240.06嚴(c)實驗題目五curve fittingleast mean square error姓名:xxxxxx學號：xxxxxxxxxxxx時間：2015年5月28 口【題目】1 the temperatu

34、re cycle in a suburb of los angeles on november 8 is given in theaccompanying table. there are 24 data points.time. pjn.degreestime«ajn.degrees1661582662*583653584644585635576636577627578618589609601060106411591167midnight58noon68(a) follow the procedure outlined in example5.5 (use the fmins co

35、mmand) to find the least-squares curve of the form f (x) = a cos(bx)+c sin(z)x)+ e for the given set of data.(b) determine e2(f).(c) plot the data and the least-squares curve from part (a) on the same coordinate system.【實驗分析】(a) in this problem the hour 1-midnight p.m. can be map to 1 12, 1 noon a.m

36、. can be map to 1324. so hour and temperature can be used to fit the curve to get the curve with minimum error, which is define as:nnerror = 工(f(xk) 一 y/c)2 =g4cos(b$) + csin(dxj + e - yk)2 i=li=lthe best a, b, c, d, e can be found when:derror derrorderror derrorderrordadbdcddde(b) e2 (/) is defined

37、 as:nk=lin our case,(n-acos(ibxk) + csin(dxk) + e - yk2 k=l【實驗過程和程序】(a) the function of the error is:function z=e_5(u)x=l:24:y=66» 66, 65, 64, 63, 63, 62,61, 60, 60, 59. 58, 5& 58, 38, 58. 57, 57, 57, 58, 60. 64, 67, 68: a-u(l);b-u(2)：cu：d-u(4)：e-u(5)：z»0：for i=l:length(x)z=z+ (a*cos (

38、b*x(i)-k：tsin(d*x(i)+e-y(i) "2;endusing command fminsearch< e-? »-11 1 1 o0- a b c d e can be solved.(b)by the function solved by part(a), we can calculate /(%)fdr each sample point andthen calculate e2 (/):x=l:24;y-：66, 66, 65, 64, 63, 63, 62, 61,60, 60, 59, 5& 58, 5& 5& 58

39、, 57, 57, 57, 58, 60, 64, 67, 68：:para= fminsearchc e_5', 1 0.1 1 1 60);e2f=0;for i=l:length(x)e2f«e2f+(para(l) »cos (para(2) *x (1)+para(3) »sm(para(4) *x (1) )+para(5)-y(i) 2/length(x): ende2f«e2f*0.5(c) code to plot the line:x=l:24：y=66, 66, 65, 64, 63, 63, 62, 61. 60, 60,

40、 59, 58, 3& 58, 58, 58, 57r 57, 57, 58, 60, 64, 67, 68: para= fninsearchc e_5 , 1 0.1 1 1 60：);x_new»1:0. 01:24:y_new«para(l)*cos(para(2)*x_new)+p ar a(3)»sin(para(4)*x_new)+para(5);plot (x, y, ' r*'»x.new, y_new):【實驗結(jié)果】(a) the result for parameter ab c d e is:ans =4.

41、86720.2213 0.31791.2896 61.9214so the least-squares curve is :/(%) = 48672cos(0.2213心)一 0.3179sm(1.2896xk) + 61.9214(b) error for f(x) is:e2f 1. 4349(c) graph for least-squares curve and sample points is:notice that the result of function fminesarch is sensitive about the input seeds. firstly the se

42、ed is 1 1 1 1 1, and get the result 0.1055 -2.1100 0.9490 -4.8698 61.0712. in second step, seed is 0.1 2 1 -5 61 and result is 0.0932 2114 175 5127 61.0424. by 4 iterations, i find a good fit of the data and get the seed 1 0.1 1 1 60.實驗題目六numerical integrationautomatically select the integration ste

43、p of the trapezoidal method姓名：xxxxxx學號：xxxxxxxxxxxx時間：2015年5月28 口【題目】(i) approximate each integralusing the composite(a)(b)/l(i +x2)-v£r (2 + sin(2/r)dztrapezoidal rule with m= 10.j-ijo(c)(d)(ii) approximate each integralusing the compositejo.25 vxjosimpson rule with m= 5.【實驗分析】(e)(f)f 2xcas(j)

44、dx1 sin(2x)exdxjojothe formula of each function is given in this problem, so we just need to usecomposite trapezoidal rule and composite simpson rule to solve them.【實驗過程和程序】function of the six equations:function result=fl (x)function result=f2(x) function result=f3(x) result=(l+xa2)a(-1) ; result=2+

45、sin(2*x 0.5): result=l/x*0.5： function result=f4(x) function result=fo(x) function result=f6(x) result=xa 2*exp(-x); result=2*x*cos(x); result=sin(2*x)*exp(-x);code for calculation:trap, trap trap trap trap trapsuml=traprl(f 1,-1,sum2=traprl(f2,0,sum3=traprl (f 3, 0. 25, sum4=traprl(f4, sumo=traprl(f5, sum6=traprl(f6,0,0,1,10)1,10)4,10)4.10)2.10)0, pi, 10)simp simp simp simp simp simpsuml=simprl(f1,-1sum2=simprl(f2,0sum3=simprl(f3, 0. 25suml=simprl