混合週期因子的生長曲線建模 - 政大學術集成

全文

(1)國⽴政治⼤學應⽤數學系碩⼠學位論⽂. 立. 政治大. 混合週期因⼦的⽣⾧曲線建模 ‧ 國. 學 ‧. The Modeling of Logistic Curve by Mixing the Periodic Factor n. er. io. sit. y. Nat. al. Ch. engchi. i n U. v. 碩⼠班學⽣：林昱撰指導教授：曾正男博⼠中華民國 106 年 5 ⽉ 11 ⽇.

(2) 國⽴政治⼤學應⽤數學系林昱君所撰之碩⼠學位論⽂混合週期因⼦的⽣⾧曲線建模 The Modeling of Logistic Curve by Mixing the Periodic Factor 政治大立. ‧. ‧ 國. 學. 業經本委員會審議通過論⽂考試委員會委員：. n. er. io. sit. y. Nat. al. Ch. engchi. i n U. v. 指導教授：系主任：. 中華民國 106 年 5 ⽉ 11 ⽇.

(3) 謝耶! 終於把論⽂打完囉!! 真是太開勳了!! 這些⽇⼦要感謝曾正男⽼師對我的協助與指導，讓我能夠在寫這篇論⽂時，能夠⼗分得⼼應⼿。同時也要感謝幫我⼝試的薛⽼師和李⽼師，對我的論⽂提出許多修改的建議。. 政治大. 在碩班的這幾年，真是過得多采多姿。⾮常幸運的在這段期間遇到了我的⼥朋友⾠兒，⼀起在啦啦隊裡同⽢共苦，參加了⼤⼤⼩⼩的⽐賽跟表演。在系上，也交到了⼀群⽩痴朋友阿宗、治安、沛倫跟⼩潤，⼀起度過了實變的激戰還有各種夜晚的促膝⾧談。除此之外，這段期間也很開⼼加⼊了政⼤啦啦隊，能跟隊上的⼤家⼀起在場上揮灑青春的汗⽔，應該是我碩⼠班⽣涯中美好的⼀段回憶。最後，也很⾼興能跟屁孩⼀起踏⼊重訓的這個領域，希望未來能⼀起練成超強巨巨。. 立. ‧. ‧ 國. 學. sit. y. Nat. n. al. er. io. 能完成碩⼠學位算是⼈⽣中重要的⼀個⾥程碑。接下來的⽇⼦希望能將所學做出對社會有貢獻的事。. Ch. engchi. i n U. v. 林昱謹誌于國⽴政治⼤學應⽤數學所中華民國⼀百零六年七⽉. ii.

(4) 中⽂. 要. 在這篇論⽂中，我們引⽤“海岸綠提-⽔筆仔”的資料。我們的⽬的是要找到⼀個合適的函數去擬合這些資料，以及找到⼀個合適的微分⽅程式，使得該⽅程式的解是能夠擬合這些資料的。這樣的函數以及微分⽅程將能幫助我們更了解這筆資料的性質，並對預測資料未來⾛向更有幫助。. 政治大我們⾸先藉由⽣⾧曲線來建構我們的數學函數，並對這個數學函數的模立 ‧. ‧ 國. 學. 型加上週期項來改良。藉由 Matlab curve fittig tool，我們找到了這個函數的⼀組參數來擬合原始資料。最後得到的結果如我們的預期，加了週期函數做出來的建模較原本⽣⾧曲線的建模更為貼近原始資料。. n. al. er. io. sit. y. Nat. 在尋找微分⽅程系統的建模上，我們參考了⽂獻⽅法，並加以改良。然⽽這個新的微分⽅程式在加上週期項來改良之後卻沒辦法找到解析解，所以我們利⽤基因演算的⽅法來去尋找適合這個微分⽅程系統的參數，並搭配 Heun’s method，RK2 method 和 RK4 method 求出⼀些數值解。最後得到了⼀個⽐原參考⽂章更好的結果。. Ch. engchi. i n U. v. 關鍵字：曲線擬合、數學模型、⽣⾧曲線、基因演算法. iii.

(5) Abstract In this paper, we focus on the data from the website ‘Seacoast Green BankKandelia’. There are two things we want to do for these data. First, we want to. 政治大 equation such that its solution 立 fits these data well. By exploring the function and find a function which graph fits these data. Second, we want to find a differential. ‧ 國. 學. the differential equation, we can understand more properties of these data.. ‧. We first build our mathematical function from Logistic curve and improve it by. y. Nat. adding a periodic factor. By using Matlab curve fitting tools, we find parameters. io. sit. of this function which fit these data well. The final result is the same as our expec-. n. al. er. tation. The model by adding a periodic factor fits the data better than the model of Logistic function.. Ch. engchi. i n U. v. To look for a differential equation, we follow the method in [7] and improve it. However, there is no analytical solution after adding a periodic factor into the model. Thus we use the method of a genetic algorithm to find suitable parameters of this differential equation. Moreover, we find the numerical solution by using Heun’s method, RK2 method and RK4 method. Finally, we get a better result than one in Ren-fa, Chen’s paper.. Keywords: Curve fitting, Mathematical model, Logistic curve, Genetic algorithm. iv.

(6) Contents ⼝試委員會審定. i. 謝. 立. 要. sit. n. al. er. io List of Tables. iv v. y. Nat. List of Figures. ‧. Contents. 學. Abstract. ii iii. ‧ 國. 中⽂. 政治大. Ch. engchi U. v ni. vii ix. 1. Introductin. 2. Source of Data. 4. 2.1. Sources . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 4. 2.2. Time and Height . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 7. 3. 1. Construct Model Function. 9. 3.1. Fit by Logistic Curve . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 9. 3.2. Model 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 10. 3.3. Model 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 11. 3.4. Genetic Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 13. v.

(7) Construct Differential Equation Model. 15. 4.1. Survey from Ren-fa’s Paper . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 15. 4.2. First Order Differential of the Data . . . . . . . . . . . . . . . . . . . . . . . .. 20. 4.3. Model 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 22. 4.3.1. Using Regression . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 23. 4.3.2. Using Finite Difference Method . . . . . . . . . . . . . . . . . . . . .. 25. 4.3.3. Using Analytical Solution . . . . . . . . . . . . . . . . . . . . . . . .. 29. Model 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 32. 4.4.1. Using Finite Difference Method . . . . . . . . . . . . . . . . . . . . .. 32. 4.4.2. Using Genetic Algorithm . . . . . . . . . . . . . . . . . . . . . . . . .. 35. 4.4.3. Using Matlab Curve Fitting Tool . . . . . . . . . . . . . . . . . . . . .. 4.4. Conclusion. 學. 5. 立. 政治大. ‧. ‧ 國. 4. A Code use in paper. sit. y. Nat. A.1 Code 01 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. io. er. A.2 Code 02 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.3 Code 03 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. n. al. Ch. n U i engch. iv. 39 42 44 44 55 65. A.4 Code 04 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 76. A.5 Code 05 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 90. Bibliography. 102. vi.

(8) List of Figures 2.1. Some pictures of process of growth in different month . . . . . . . . . . . . . .. 7. 2.2. Graph of time and height . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 7. 3.1. 政治大 Graph of data and equation 3.2.1 . . . . . . . . . . . . . . . . . . . . . . . . . 立. 11. Graph of data and (3.3.1) . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 12. 3.3. Predict of Kandelia’s height by equation 3.3.1 . . . . . . . . . . . . . . . . . .. 13. 4.1. Graph of Table 2.3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 17. 4.2. Graph of Table 4.2 and equation4.1.3 . . . . . . . . . . . . . . . . . . . . . . .. 18. 4.3. Graph of E(t) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 4.4. Numerical result of equation 4.1.5 . . . . . . . . . . . . . . . . . . . . . . . .. 20. 4.5. . . . . . . . . . . . . . .. 24. ‧. ‧ 國. 學. 3.2. er. io. sit. y. Nat. al. 19. Best result . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 25. 4.7. Best result . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 29. 4.8. Graph of data and (4.3.1) . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 30. 4.9. Graph of data and (4.3.5) . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 31. 4.10 Best result . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 35. 4.11 Fit the data by genetic algorithm . . . . . . . . . . . . . . . . . . . . . . . . .. 36. 4.12 Fit the data by genetic algorithm with parameters adjusted . . . . . . . . . . .. 37. 4.13 Best result . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 38. 4.14 Predict of Kandelia’s height by equation 4.4.1 . . . . . . . . . . . . . . . . . .. 39. 4.15 Graph of data and equation(4.4.1) . . . . . . . . . . . . . . . . . . . . . . . .. 40. n. 4.6. v i n Fit the data by regression . . . . . C .h. . . . . . . . . U engchi. vii.

(9) 4.16 Best result . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 立. 政治大. ‧. ‧ 國. 學. n. er. io. sit. y. Nat. al. Ch. engchi. viii. i n U. v. 41.

(10) List of Tables Germination under different environments . . . . . . . . . . . . . . . . . . . .. 5. 2.2. Process of growth of Kandelia all the year round. 6. 2.3. . . . . . . . . . . . . . . . . 治政 Relation of time and height . . . . . . . . . .大 . . . . . . . . . . . . . . . . . . 立. 4.1. Relation between velocity of growth and height of hypocotyl of Kandelia . . .. 16. 4.2. Relation between acceleration of growth and height of hypocotyl of Kandelia .. 17. 4.3. Relation between velocity of growth and height of hypocotyl of Kandelia . . .. 21. 5.1. Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 43. ‧. ‧ 國. 學. 2.1. n. er. io. sit. y. Nat. al. Ch. engchi. ix. i n U. v. 8.

(11) Chapter 1 Introductin 政治大 In this thesis, we focus on the data from the website ‘Seacoast Green Bank-Kandelia’ [1]. 立. ‧ 國. 學. There are two things we want to do with these data. First, we want to find a function which graph fits these data. Second, we want to find a differential equation such that its solution fits. ‧. these data well. By exploring the function and the differential equation, we can understand more properties of these data.. y. Nat. io. sit. We begin our paper in chapter 2 by referring to the website [1]. In this website, we observe. n. al. er. flower of Kandelia. Also, we record the process of its anthesis and seed-bearing. Furthermore,. i n U. v. we are interested in the growth of hypocotyl. Thus, we sort out height of Kandelia in different. Ch. engchi. month. We then notice that these data perform quite similar to logistic curve. Therefore, our modeling starts by using it. Here is a brief introduction to the logistic curve according to the website [4]. Further detail will be shown in chapter 2.. A logistic function or logistic curve is a common ”S”shape (sigmoid) curve, using the equation below: Q(t) =. B 1 + Ae− Bkt. Here, e is the natural logarithm (also known as Euler’s number). A , B and k are constants.. Next, in chapter 3, we build our mathematical function based on the logistic curve. By using. 1.

(12) the Matlab curve fitting tool, we find suitable parameters. We want the function with these parameters fit original data well. The model is below: Let Q(t) be the height of Kandelia, A , B and k are the parameters.. Q(t) =. B 1 + Ae− Bkt. We then fit our data to this model. After that, we notice that the curve seems to be different from the original data by adding a periodic factor. Thus we adjust our model as follows:. 政治大. Again, let Q(t) be the height of Kandelia. A , B, k, a, b, and c are the parameters we want. 立. Q(t) =. 學. ‧ 國. to find.. B + asin(bt + c) 1 + Ae− Bkt. ‧. The final result is the same as our expectation. By adding a periodic factor, the model fits. Nat. sit. y. the data better than the model of the logistic function. However, the parameters suitable for. er. io. this model were found by using Matlab curve fitting tool. The problem is, we are not familiar. al. v i n C h way of gettingUthese parameters. Genetic Algorithm [6]. This is a clearer engchi n. with the operation inside the Matlab. Thus, at the end of Chapter 3, we give an introduction to. Next in chapter 4, to look for a differential equation model, we follow R-f’s paper [7]. We then improve it by estimating the velocity of growth of the Kandelia in a better way. That is, instead of using the forward difference method [5] R-f’s paper does, we use the central difference method. This will make estimation much better. Note that the main difference between the models in chapter 3 and chapter 4 is that, in chapter 3 we find a function to fit the data directly. On the other hand, in chapter 4 we alternate the form of the logistic curve into a first order differential equation. We want to find parameters which fit the differential equation. Thus we construct the following model: Q ′ ( t ) = p1 Q2 ( t ) + p2 Q ( t ). 2.

(13) Here, Q′ (t) is the growth rate of Kadelia. Q(t) is the height of Kandelia. Also, p1 and p2 are the parameters we hope to find.. Under this model, we use the relationship between height and velocity. With regression, we can find the suitable parameters. Next, by using Heun’s method, the RK2 method and the RK4 method, we then find the numerical solution of this model. Furthermore, since this equation is the form of Bernoulli equation [2], we can find its analytical solution. We then compare this solution to the numerical result.. 政治大. Similar to Chapter 3, we notice that the curve seems to be different from the original data. 立. after adding a periodic factor. Thus we adjust our model as following:. ‧ 國. 學. Q′ (t) = p1 Q2 (t) + p2 Q(t) + asin(bQ(t) + c). ‧. Q′ (t) is the growth rate of Kadelia. Q(t) is the height of Kandelia. Here, a, b, c, p1 and p2. y. Nat. n. er. io. al. sit. are the parameters we hope to find.. i n U. v. However, there is no analytical solution to the model. Thus we use genetic algorithm [6] to. Ch. engchi. find suitable parameters for this differential equation. Moreover, we find the numerical solution by using Heun’s method, the RK2 method and the RK4 method. Finally, we get a better result than the survey paper [7]. In chapter 5, we give a review to all the methods we used in this paper. Also, we draw a table to show all the results.. 3.

(14) Chapter 2 Source of Data 政治大 In this chapter, we will talk about the data we use in this paper. 立. In Section 2.1, we will. ‧ 國. 學. give the website [1] an introduction. It is a website recording behaviors of Kandelia. Next, in Section 2.2, we will sort out the data of Kandelia from time to time.. ‧ y. sit. Sources. Nat. 2.1. n. al. er. io. We focus on two observations of Kadelia in the website [1]. First, we observe hypocotyl of. i n U. v. Kandelia in four different environment. We then recored its phenomenon of germination in six weeks.. Ch. engchi. Second, we observs flower of Kandelia. Also, we record the process of its anthesis and seedbearing. Furthermore, we are interested in the growth of hypocotyl. Detailed observations are shown in Table 2.1 and Table 2.2. We also put some pictures of process of growth in different months. From left to right namely: “small bud in March”, “most of flower blossom in July”, “seed in September”, “fruit grow viviparous seedlings in October”, and “viviparous seedlings grow about twenty centimeters February next year”. The pictures are shown in Figure 2.1.. 4.

(15) 5. Viviparous seedlings grow small roots at the end.. Viviparous seedlings grow small roots at the end.. Insert into water. Sink into water. Soil of Intertidal zone and general water. Viviparous seedlings grow about 0.4 to 0.6 centimeters of root.. Viviparous seedlings grow about 0.5 to 0.6 centimeters of root.. Soil of Intertidal zone and water of Intertidal zone. General soil and general water. First week Viviparous seedlings grow about 0.5 to 0.7 centimeters of root.. er. sit. Second week Viviparous seedlings about 1.5 to 2 centimeters of root,end of bud begins to grow. Viviparous seedlings about 1.5 to 2 centimeters of root,end of bud begins to grow. Viviparous seedlings about 1.5 to 2 centimeters of root,end of bud begins to grow. Viviparous seedlings grow small roots at the end, bud end almost no growth. Viviparous seedlings grow about 0.3 centimeters of small roots at the end.. n. Ch. engchi. There are eight viviparous seedlings start to grow leaves.. There are eight viviparous seedlings to sprout of phenomena.. i n U. Viviparous seedlings grow about 1 centimeters of small roots at the end.. end almost no growth.. Viviparous seedlings grow. ysmall roots at the end, bud. Viviparous seedlings grow about 1.5 to 2 centimeters of small roots at the end but a few number of roots.. Viviparous seedlings grow small roots at the end, bud end almost no growth.. Viviparous seedlings grow about 1.5 to 2 centimeters of small roots at the end but a few number of roots.. Viviparous seedlings grow small roots at the end, bud end begins to grow a little.. There are nine viviparous seedlings grow a pair of leaves.. There are five viviparous seedlings grow a pair of leaves.. There are four viviparous seedlings start to grow leaves.. There are four viviparous seedlings to sprout of phenomena.. ‧ 國. Sixth week There are four viviparous seedlings grow a pair of leaves.. Fourth week There are two viviparous seedlings start to grow leaves.. Third week There are three viviparous seedlings to sprout of phenomena.. ‧. io. al. 學. Nat. Table 2.1: Germination under different environments. 立政治大. v.

(16) Table 2.2: Process of growth of Kandelia all the year round Month Last year of February. March. April. Observed Results Hypocotyls of Kandelia grow about eighteen to twenty-five centimeters, at the end of February a little part of Kandelia began to fall. Kandelia grow new small bud like small rod. Hypocotyls of Kandelia grow about eighteen to twenty-five centimeters, at the end of March most Kandelia began to fall and start defoliating. Small buds of Kandelia grow in a shape resembling a matchstick. 1. Previously falling Kandelia viviparous seedlings sprout almost together at the middle of April, we found Hypocotyl sprout has nothing to do with length,the main factor is to be fixed hypocotyl about 5-10 cm deep in sediment, hypocotyl not to be washed away or moved with tide and river. 2. Hypocotyls of Kandelia grow about twenty to twenty seven centimeters. Kandelia almost fall at the end of April. Kandelia and small rod small bud begin to grow into the shape of a small fire wood stick. Kandelia viviparous seedlings grow root in two weeks and grow new leaves small rod small bud about 1 cm in two months. One part of Kandelia blossom and the other part of Kandelia is also bud. All Kandelia blossom, star-shaped flowers with a fragrance and attract a large number of foraging. Brown seeds grow. Small fruits grow one centimeter in width and two centimeter in length, at the end of September all the seeds will develop a small hypocotyl. Hypocotyl of Kandelia grow about one to eight centimeters. Hypocotyl of Kandelia grow about eight to fifteen centimeters. Hypocotyl of Kandelia grow about ten to twenty centimeters. Hypocotyl of Kandelia grow about fifteen to twentytwo centimeters. Hypocotyl of Kandelia grow about eighteen to twentyfive centimeters.. 立. 政治大. ‧ 國. November December January February. ‧. October. Kandelia (small rod of small bud) begin to grow into the shape of a small fire wood stick.. 學. August September. y. sit. er. al. n. July. io. June. Nat. May. Note Dropping length of hypocotyl short than the average hypocotyl about four centimeters. At the end of March almost all of hypocotyls of Kandelia have fallen. Ch. engchi. 6. i n U. v. A small part of bud will blossom at the end of May. The end of June is the peak of flowering. Blossom ends in midJuly. Size as peanuts. Length of Viviparous seedlings vary. Length of Viviparous seedlings vary. Length of Viviparous seedlings vary. Length of Viviparous seedlings vary. Length of Viviparous seedlings vary. Length of Viviparous seedlings vary..

(17) Figure 2.1: Some pictures of process of growth in different month. 2.2. Time and Height. Accroding to Section 2.1, we know the growth of hypocotyl of Kandelia weekly and monthly.. 治政大 Further data is in Table 2.3 and to height of growth each months from September to February. 立 relations are shown in Figure 2.2.. We will reasonably assume the height of growth each weeks. In addition, we give restriction. ‧. ‧ 國. 學. n. er. io. sit. y. Nat. al. Ch. engchi. i n U. v. Figure 2.2: Graph of time and height. 7.

(18) Table 2.3: Relation of time and height. io. n. Ch. January second week January third week January fourth week February first week February second week February third week February fourth week March first week March second week March third week March fourth week April first week April second week April third week April fourth week. engchi. 8. y. sit. Nat. al. 0.1 0.3 0.6 0.8 1 2.5 3 4.5 5.7 7.9 9.2 11.5 12 13.2 14.2 15. ‧. ‧ 國. 立. 學. Initial September first week September second week September third week September fourth week October first week October second week October third week October fourth week November first week November second week November third week November fourth week December first week December second week December third week December fourth week. 治 Time 政Height u(t) 大 0 January first week. Period t 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16. er. Time. i n U. v. Period t 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32. Height u(t) 15.7 16.8 17.5 18.5 19.25 20.1 20.85 21.5 22.15 22.75 23.15 23.5 24 24.4 24.75 25.

(19) Chapter 3 Construct Model Function 政治大 In this chapter, we first introduce logistic curve in section 3.1. Then we construct a model 立. ‧ 國. 學. function by logistic curve in section 3.2. Also,we improve this model in section 3.3. Whether we use logistic model or the improved model, there are many unknown parameters we need to. y. Nat. io. n. al. sit. Fit by Logistic Curve. er. 3.1. ‧. find. Therefore, in section 3.4 we will introduce Genetic algorithm.. i n U. v. We begin by observing the data in Figure 2.2. Notice that these data is increasing under time.. Ch. engchi. The increasing rate grows slow from beginning. Then, it becomes faster in the middle. After that, it slow down again. The performs of these data is quite similar to logistic curve [4]. Thus we give logistic curve an introduction here: A logistic function or logistic curve is a common ”S” shape (sigmoid curve), with equation 3.2.1 Q(t) =. B 1 + Ae− Bkt. where. • A , B and k are constants • e = the natural logarithm base (also known as Euler’s number) 9.

(20) The function was named in 1844–1845 by Pierre François Verhulst, who studied it in relation to population growth. The initial stage of growth is approximately exponential; then, as saturation begins, the growth slows, and at maturity, growth stops. The logistic function finds applications in a range of fields, including artificial neural networks, biology (especially ecology), biomathematics, chemistry, demography, economics, geoscience, mathematical psychology, probability, sociology, political science, linguistics, and statistics.. 3.2. Model 1. 立. 政治大. ‧ 國. 學. We begin our modeling by using logistic curve mentioned in section 3.1. That is, we are going to use the data in Table 2.3 as the value of Q(t). We want to find parameters A, B, and k. (3.2.1). er. io. al. B 1 + Ae− Bkt. sit. Nat. Q(t) =. y. ‧. which fit the following equation:. v. n. By using Matlab curve fitting tool, we then find A , B and k as follow: • A = 25.8. Ch. engchi. i n U. • B = 24.21 • k = 0.009621 Substituted the above parameters into equation 3.2.1, we have result shown in Figure 3.1.. 10.

(21) Figure 3.1: Graph of data and equation 3.2.1. 政治大 By this method, the sum of square error (SSE) is “5.38952958”. Overall, this is not the ideal 立 result we want. Thus, a further improvement will be done in next section.. ‧ 國. 學. Model 2. ‧. 3.3. Nat. sit. y. We now observe the curve in Figure 3.1. The value of model function is higher then the data. er. io. from week 0 to week 8. Then from week 9 to week 16, its value is lower then the data. Next,. al. v i n C h seems to be different 32. The behavior of this model function e n g c h i U from the original data by adding a n. from week 17 to week 24 it goes higher again. After that, it goes lower from week 25 to week. periodic factor. Thus, in order to decrease the SSE, we now consider equation 3.3.1 as our new model: Q(t) =. B + asin(bt + c) 1 + Ae− Bkt. (3.3.1). Similar to section 3.2, we use the data in Table 2.3 as the value of Q(t) and find parameters a, b, c, A, B, and k which fit the following equation well. By using curve fitting tool again, we have the parameters as below: • A = 29.83 • B = 23.79 • k = 0.01033 11.

(22) • a = 1.233 • b = 0.3468 • c = 8.948 We substitute the parameters above into equation 3.3.1. The result is shown in Figure 3.2.. 立. 政治大. ‧. ‧ 國. 學 Figure 3.2: Graph of data and (3.3.1). sit. y. Nat. io. er. By this method, the SSE is “4.92591951”. That is a pretty good result. Also, compare this result to the one in section 3.2, we have:. n. al. i n U. v. C h4.92591951 i 1 − e n g c h≈ 0.08602 5.38952958. which means we reduced about 8.6% of sum of square error. The model we built here seems to fit pretty well. However, its actually over fit our data. To be more specific, we take a closer look at equation 3.3.1 again. As the time keeps going, we plot the graph of prediction:. 12.

(23) 政治大. Figure 3.3: Predict of Kandelia’s height by equation 3.3.1. 立. ‧ 國. 學. In Figure 3.3. We see the curve remains oscillate after t = 30. Since this is a data from height of Kandelia, this kind of behavior is not reasonable. In the real world, the height is more. ‧. likely to keep growing. Also, it will end up approaching to a stable value. Therefore, we will. sit er. io. 3.4. y. Nat. try improving this part in model later.. Genetic Algorithm a. n. iv l C n h curve i Ufind the parameters of models in sece n gfitting We are curious about how Matlab c h tools. tion 3.2 and section 3.3. With some survey, we find it works similar to the process of genetic algorithm [6]. Thus, we will give an introduction to genetic algorithm here. This algorithm will be used frequently in the follow up chapter.. In computer science and operations research, a genetic algorithm (GA) is a metaheuristic inspired by the process of natural selection that belongs to the larger class of evolutionary algorithms (EA). Genetic algorithms are commonly used to generate high-quality solutions to optimization and search problems by relying on bio-inspired operators such as mutation, crossover and selection.. 13.

(24) Now given data ( xi , yi )in=0 . We construct a model function y = f ( a1 , a2 , ..., am , x ), where a1 , ..., am are the parameters and x ∈ [ x0 , xn ]. Our goal is to find a set of a1 , ..., am , such that the function with these parameters has a relative small sum of square error with respect to the data.. (0). (0). Give an initial set a1 , ..., am . With these parameters, we can compute the sum of square (0). error SSE(0) . Let ai = ai , i = 1...m, and SSE = SSE(0). Next, we give each ai ,i = 1...m, some perturbation and get a set of new parameters (1). (1). a1 , ..., am . Again, we can compute another sum of square error SSE(1) . Compare SSE and SSE(1) :. 立. (1). In this case, we replace each ai by ai (1). and replace SSE by SSE(1) . Denote by:. for i = 1...m. ‧. ai = ai. 學. ‧ 國. • case1:SSE > SSE(1). 政治大. SSE = SSE(1). y. Nat. er. io. Here we do nothing. sit. • case2:SSE ≤ SSE(1). al. n. v i ( 2 ) (2) n Again we give our new a1 , ...,C amhperturbation and U get a1 , ..., am . As before, we find the i engch (2). corresponding sum of square error SSE. . Similar to the former two cases, we renew ai for. i = 1...m and SSE.. Finally, repeating above steps, we can find a set of parameters a1 , ..., am which corresponding SSE is relative small.. 14.

(25) Chapter 4 Construct Differential Equation Model 政治大 In this chapter, we want to look for a differential equation as our model. A reference of R-f’s 立. ‧ 國. 學. paper is shown in section 4.1. Then, in section 4.2, we will improve his method of estimating Kandelia’s growth velocity. Next, in section 4.3, we construct a model by changing Logistic. ‧. function into form of first order differential equation. After that, in section 4.4, we improve this model by adding a periodic factor.. Survey from Ren-fa’s Paper. n. al. Ch. engchi. er. io. sit. y. Nat. 4.1. i n U. v. We want to find a differential equation model. To do that, we follow R-f’s paper [7]. Here is the idea. First, we find growth velocity of Kandelia from the data of its height. Next, we use Kandelia’s growth velocity to find acceleration rate of it. Thus, we have the acceleration rate of Kandelia in different height. In other word, we can plot a graph with Kandelia’s height as x axis and Kandelia’s acceleration rate as y axis. Our goal is to find an equation that fits this graph. We give a summary of R-f’s work:. In this paper, we first calculate the first order differential value of the data in Table 2.3 by using forward difference approximation. These data will represent the growth velocity of Kandelia in different height:. 15.

(26) h ′ ( ti ) ≈. h ( t i +1 ) − h ( t i ) t i +1 − t i. (4.1.1). Here, h′ (ti ) is the approximation of Kandelia’s growth velocity at time ti , and h(ti ) is the height of Kandelia at time ti . The result is shown in Table 4.1: Table 4.1: Relation between velocity of growth and height of hypocotyl of Kandelia Velocity of growth 0.1 0.2 0.3 0.2 0.2 1.5 0.5 1.5 1.2 2.2 1.3 2.3 0.5 1.2 1.0 0.8 0.7. 政治大. n. Ch. engchi. Velocity of growth 1.1 0.7 1.0 0.75 0.85 0.75 0.65 0.55 0.6 0.15 0.4 0.35 0.4 0.35 0.25. y. sit. er. io. al. Height 15.7 16.8 17.5 18.5 19.25 20.1 20.85 21.5 22.15 22.75 23.15 23.5 24 24.4 24.75 25. ‧. ‧ 國. 立. Period 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32. 學. Height 0 0.1 0.3 0.6 0.8 1 2.5 3 4.5 5.7 7.9 9.2 11.5 12 13.2 14.2 15. Nat. Period 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16. i n U. v. Next, we use forward difference approximation on Table 4.1. Thus, we can calculate the second order differential value: h′′ (ti ) ≈. h ′ ( t i +1 ) − h ′ ( t i ) t i +1 − t i. (4.1.2). Here, h′′ (ti ) is the growth acceleration of Kandelia at time ti , and h′ (ti ) is the growth velocity of Kandelia at time ti . The result is shown in Table 4.2.. 16.

(27) Table 4.2: Relation between acceleration of growth and height of hypocotyl of Kandelia Acceleration of growth 0.1 0.1 -0.1 0 1.3 -1 1 -0.3 1 -0.9 1 -1.8 0.7 -0.2 -0.2 -0.1 0.4. 立. Period 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32. Height 15.7 16.8 17.5 18.5 19.25 20.1 20.85 21.5 22.15 22.75 23.15 23.5 24 24.4 24.75 25. 政治大. Acceleration of growth -0.4 0.3 -0.25 0.1 -0.1 -0.1 0 -0.05 -0.2 -0.05 0.15 -0.1 -0.05 -0.1. ‧. ‧ 國. Height 0 0.1 0.3 0.6 0.8 1 2.5 3 4.5 5.7 7.9 9.2 11.5 12 13.2 14.2 15. 學. Period 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16. sit. y. Nat. Once we have the data in Table 4.2 and Table 2.3, we get the relation between second order. io. er. differential data and original data. We then have a graph with Kandelia’s height as x axis and its acceleration rate as y axis. It is shown as in Figure 4.1. We want to find a curve to fit it.. n. al. Ch. engchi. i n U. v. Figure 4.1: Graph of Table 2.3 17.

(28) In R-f’s paper, we conclude that equation 4.1.3 is the best mathematical model fitting Figure 4.1: h′′ (t) = a1 sin(b1 h(t) + c1 ) + a2 sin(b2 h(t) + c2 ) + a3 sin(b3 h(t) + c3 ). (4.1.3). It is a second order differential equation with sum of sin functions. We want to find suitable parameters: ai , bi , and ci (i = 1, 2, 3). By using matlab curve fitting tool, we find result as below: • a1 = 0.4573 • b1 = 3.497. 立. ‧ 國. 學. • c1 = −1.402. 政治大. • a2 = 0.3126. y. sit er. al. n. • b3 = 2.424. io. • a3 = 0.3544. Nat. • c2 = 1.053. ‧. • b2 = 2.662. Ch. engchi. i n U. v. • c3 = 2.24 We then substitute above parameters into equation 4.1.3. The result is shown in Figure 4.2.. Figure 4.2: Graph of Table 4.2 and equation4.1.3 18.

(29) Next, we try to find an analytical solution of equation 4.1.3. The result is as follow: ( ) 3 ′′ If h(t) is the solution of h (t) = ∑ ai sin bi h(t) + ci and t0 = 0, h(t0 ) > 0, then we i =1. have : 1. ( ) 3 1 ′ 2 a h (t) + ∑ i cos bi h(t) + ci = E(t). 2 b i =1 i. (4.1.4). 2. v u ( ( )) 3 u ai ′ t h (t) = 2 E − ∑ cos bi h(t) + ci , b i =1 i. 立. 政治大. (4.1.5). ‧ 國. 學. Here, E = E(t) is a constant. Also, ai , bi and ci are constants. We find out R-f’s paper end up here. However, we are curious about the solution of equation 4.1.5. Therefore, we use. Nat. sit. y. ‧. numerical method to find a solution of it. Here is how it works:. al. er. io. Since it has been proved that E(t) is a constant. Our first step is to find the value of it. We. v. n. use data in Table 4.1. Also, we substitute Kandelia’s velocity into h′ (t) and its corresponding. Ch. engchi. i n U. height into h(t). Thus, we have a value of E(t) in equation 4.1.4 as follow:. Figure 4.3: Graph of E(t) 19.

(30) We take a closer look at Figure 4.3. We see that E(t) is stable at t > 12. Thus, we take average of them and we have E = ”0.266860291408”. Next, we use RK4 method to solve equation 4.1.5. The result is shown below:. 立. 政治大. ‧ 國. 學. Figure 4.4: Numerical result of equation 4.1.5. ‧ sit. y. Nat. The blue line in Figure 4.4 is the original data of Kandelia’s height, and the red line is our. io. er. fitting curve. We compute the sum of square error and find SSE = “74.3052894517”.. al. n. v i n Ch 4.2 First Order Differential the Data i U e n gofc h. Recall that in R-f’s paper, we try to find growth velocity of Kendelia by the height of it. The method we used there is forward difference approximation [5]. That is, Given function f ∈ C1 : R → R, x ∈ R and t be a constant. Then, we can approximate the differential value of f at point x with the step size t as follow:. D f ( x, t) =. f ( x + t) − f ( x ) t. Notice that, if f ′ ( x ) is the real value. Then, this approximation will have first order big O error. That is, we have: D f ( x, t) = f ′ ( x ) + O(t). We want to decrease the error. To do that, we generate the data by central difference approx20.

(31) imation [5]. Again, given function f ∈ C1 : R → R, x ∈ R and t be a constant. Then we can approximate the differential value of f at point x with the step size t as follow:. D f ( x, t) =. f ( x + t) − f ( x − t) 2t. Now, if f ′ ( x ) is the real value. Then this approximation will have second order big O error. That is, we have: D f ( x, t) = f ′ ( x ) + O2 (t). This is better than forward difference approximation. Now, we let h(ti ) be Kandelia’s height at the time ti . We will use central difference approximation to generate Kandelia’s growth velocity:. 立. ′. 政 h(t 治) − h(大 t ). h ( ti ) =. i −1. i +1. (4.2.1). t i +1 − t i −1. ‧ 國. 學. Here, h′ (ti ) is Kandelia’s growth velocity we approximate at the time ti . Result is shown in Table 4.3. We will use these data in follow up sections to construct our mathematical models.. ‧. n. al. Ch. 0.1 5 0.25 0.25 0.2 0.85 1.0 1.0 1.35 1.7 1.75 1.8 1.4 0.85 1.1 0.9 0.75. engchi. 21. y. Period 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32. Height 15.7 16.8 17.5 18.5 19.25 20.1 20.85 21.5 22.15 22.75 23.15 23.5 24 24.4 24.75 25. Velocity of growth 0.9 0.9 0.85 0.875 0.8 0.8 0.7 0.65 0.625 0.5 0.375 0.425 0.45 0.375 0.3. sit. Velocity of growth. er. Height 0 0.1 0.3 0.6 0.8 1 2.5 3 4.5 5.7 7.9 9.2 11.5 12 13.2 14.2 15. io. Period 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16. Nat. Table 4.3: Relation between velocity of growth and height of hypocotyl of Kandelia. i n U. v.

(32) 4.3. Model 3. We will construct a differential equation model here. Different from what R-f did in his paper, we will do the curve fitting by using the data in Table 4.3 and Table 2.3. That is, we plot a graph with Kandelia’s height as x axis and Kandelia’s growth velocity as y axis. The goal is to find a model that fits this graph well. Here is how we do it. Recall the Logistic curve in equation 3.2.1: Q(t) =. B 1 + Ae− Bkt. Here, Q(t) will be the function of height of Kandelia, A, B and k are parameters. Next, we. 政治大. differentiate both sides of this equation. The left hand side will be kept as the growth velocity. 立. 學. AB2 ke− Bkt (1 + Ae− Bkt )2 kB − k) = Q2 (t)( Q(t). Q′ (t) =. y. ‧. Nat. (4.3.1). io. er. = BkQ(t) − kQ2 (t). sit. ‧ 國. of Kandelia. On the right hand side, we will try to substitute it with Q(t). It works as below:. al. n. Since B and k are both parameters. Thus, we let p1 = Bk and p2 = −k. Our goal is to find. i n C 4.3 the best. p1 and p2 which make equation 4.3.2hfits e nTable gchi U. v. Q ′ ( t ) = p1 Q2 ( t ) + p2 Q ( t ). (4.3.2). To do this we will use method by regression in subsection 4.3.1. Next, in subsection 4.3.2, we will use the finite difference method. Also, we will improve this method by method of genetic algorithm. At last, in subsection 4.3.3, we will use an analytical approach to find the solution. This will give us a standard to compare with.. 22.

(33) 4.3.1 Using Regression In this subsection, we use regression [3] as the key method. This is a method to find the suitable parameters of a target function. Here, we use equation 4.3.2 as our target function: Q ′ ( t ) = p1 Q2 ( t ) + p2 Q ( t ). Our goal is to find suitable p1 and p2 . Notice that this is a first order differential equation. So, we need the data of Kandelia’s height versus Kendelia’s velocity. We can get this data from Table 4.3. Here is how we do it. We will start from our target function:. 政治大. ( ) Given data Q′ (ti ), Q(ti ) , i = 0, · · · , n, we have:. 立. ‧ 國. 學. Q ′ ( t i ) = p1 Q2 ( t i ) + p2 Q ( t i ). ‧. Nat. al. . . er. io. . sit. ti . To find p1 and p2 , we will construct a linear system from it:. y. Here, Q′ (ti ) is Kandelia’s growth velocity at time ti and Q(ti ) is Kandelia’s height as time. n. Q2 ( t. Q′ (t. v. . 0 ) Q ( t0 )    0 )        Q2 ( t1 ) Q ( t1 )  p1  ′       Q ( t1 )  =  .     ..  ..   ..   .  p2   .      2 Q (tn ) Q(tn ) Q(tn ). Ch. engchi. i n U. We denote the above linear system by Ax = b. Using regression, we product A T to both side of the equation. That is, A T Ax = A T b: .  Q2 ( t0 ) Q ( t0 ).   . Q2 ( t. 0). Q ( t0 ). 0). ···. Q ( t0 ). ···. Q2 ( t.      2   n )   Q ( t1 ) Q ( t1 )   p1     . ..    .  p2 Q(tn )  ..   2 Q (tn ) Q(tn ). Q2 ( t. 23.

(34) . . Q ( t0 )     2 ( t ) Q2 ( t ) · · · Q2 ( t )  Q ( t )  Q n  0 1 1   =  .   Q ( t0 ) Q ( t1 ) · · · Q ( t n )   ..    Q(tn ) . ( ) We substitute Q′ (ti ), Q(ti ) , i = 0, · · · , n by the value in Table 4.3. The solution to p1 and p2 is: • p1 = −0.00863605 • p2 = 0.21215756. 立. 政治大. We substitute above p1 and p2 into equation 4.3.2. Then, we can plot a graph with Kandelia’s. ‧ 國. 學. height as its x axis and Kandelia’s growth velocity as its y axis. The result is shown in Figure 4.5. The blue line is our original data, and the red line is our fitting curve.. ‧. n. er. io. sit. y. Nat. al. Ch. engchi. i n U. v. Figure 4.5: Fit the data by regression Using the above p1 and p2 together with the initial value in Table 2.3, we find the differential equation which fits the data of Kandelia’s velocity. We want to know whether the solution of it fit Kandelia’s height or not. Thus, using Heun’s method, RK2 method and RK4 method, we have the following result:. 24.

(35) • SSE solved by Heun’s method = 6.3916710696 • SSE solved by RK2 method = 6.39166483931 • SSE solved by RK4 method = 6.39165419709 The best result solved by this method is shown in Figure 4.6. The blue line is the data of Kandelia’s height, and the red line is our fitting curve.. 立. 政治大. ‧. ‧ 國. 學 er. io. sit. y. Nat. al. Figure 4.6: Best result. n. v i n C Notice that the best result we findhinesection 3.3 isi SSE n g c h U = “4.92591951”. It is better than the best result here: SSE = “6.3916”. Thus, we will improve our model in the follow up subsection.. 4.3.2 Using Finite Difference Method In this subsection, we will use three steps to find the suitable parameters of equation 4.3.2: Q ′ ( t i ) = p1 Q2 ( t i ) + p2 Q ( t i ). First of all, we construct the target function from equation 4.3.2. But different from section 4.3.1, we will use the finite difference method here. By this method, we can substitute Kandelia’s growth velocity into terms of Kandelia’s height. This gives us a new model with. 25.

(36) time and Kandelia’s height only. Thus, we use regression with the data in Table 2.3. We then find suitable parameters for this model. Secondly, we use genetic algorithm to find suitable initial value. That is, with the parameters above, we give the initial value some perturbation. Comparing the outcome via different initial value, we find a better result. Finally, we use genetic algorithm on the above model again. But this time, we use it not only on the initial value, but also on the parameters of the model. To be more precise, we first find an initial value and a set of parameters. Then we give the initial value some perturbation and get a new initial value. With this new initial value, we give parameters perturbation. Repeating. 政治大. this process, we will have a better outcome. The follow up is how it all works:. 立. ‧ 國. 學. From Table 2.3, we have data {( Q(ti ), ti )}in=1 , where Q(ti ) is the hight of Kandelia at time ti . Our model start from equation 4.3.2:. sit. y. ‧. Nat. Q ′ ( t i ) = p1 Q2 ( t i ) + p2 Q ( t i ). n. al. er. io. Where Q′ (t) is Kandelia’s growth velocity. Now, instead of using data in Table 4.3 on Q′ (t),. i n U. v. we use finite difference method. In other word, we let hk = tk+1 − tk−1 , (k = 1, ..., n − 1) and substituted Q′ (ti ) by. Ch. Q(ti+1 )− Q(ti−1 ) . t i +1 − t i −1. engchi. Thus, we have:. Q ( t k +1 ) − Q ( t k −1 ) = p1 Q ( t k )2 + p2 Q ( t k ) hk. ⇒ Q(tk+1 ) = Q(tk−1 ) + hk ( p1 Q(tk )2 + p2 Q(tk )). (4.3.3). The equation 4.3.3 is our target function. Notice that, with p1 and p2 be the changing vari-. 26.

(37) able, we can get the following linear system: . . h1 Q ( t1.        . . . Q ( t2 ) − Q ( t0 )         Q ( t3 ) − Q ( t1 )  h2 Q ( t2 )  p 1       =   .. ..    . .  p2      h n −1 Q ( t n −1 ) Q ( t n ) − Q ( t n −2 ). )2. h1 Q ( t1 ). h2 Q ( t2 )2 .. . h n −1 Q ( t n −1 ) 2. Denote the above linear system by Ax = b, and we will use regression to find suitable p1 and p2 . That is, we product A T to both side of the equation. The solution of the linear system A T Ax = A T b is as below:. 立. ‧ 國. ‧. . . n. al. Ch. sit. 2 2  h1 Q ( t1 ) h2 Q ( t2 ) = h1 Q ( t1 ) h2 Q ( t2 ). io. Q ( t2 ) − Q ( t0 )     2  · · · h n Q ( t n −1 )   Q ( t 3 ) − Q ( t 1 )     .  .. · · · h n Q ( t n −1 )      Q ( t n ) − Q ( t n −2 ). y. Nat. . h Q ( t1 h1 Q ( t1 )    1   2 2  · · · h n Q ( t n −1 )   h 2 Q ( t 2 ) h2 Q ( t2 )    p1     .. ..  · · · h n Q ( t n −1 )  . .   p2   2 h n Q ( t n −1 ) h n Q ( t n −1 ). 學. 2 2  h1 Q ( t1 ) h2 Q ( t2 )  h1 Q ( t1 ) h2 Q ( t2 ). . )2. er. . 政治 大. engchi. i n U. v. Using the data in Table 2.3 ,and the code in appendix A.1, we can compute p1 and p2 which fit equation 4.3.3 well. • p1 = −0.00844749 • p2 = 0.21174495 With the parameters above, we can solve equation 4.3.2 by Heun’s method, RK2 method and RK4 method. The python code is also in appendix A.1. • SSE solved by Heun’s method = 55.7242174456 • SSE solved by RK2 method = 55.66572951 27.

(38) • SSE solved by RK4 method = 55.3013480532 The result above is not very good. In order to lower the SSE, we pull in the idea of genetic algorithm. That is, we add some perturbation to the initial value. Then we can get different curves by using equation 4.3.3 within different initial value. Further detail is written in appendix A.1. Overall, we will find a curve fits Table 2.3 relatively good. With the initial value of this curve and p1 and p2 above, we have the result:. • initial value = (0.18824931457, 0.369491784562). 治政 • SSE solved by Heun’s method = 30.2074082818 大立 ‧ 國. 學. • SSE solved by RK2 method = 30.1259883441 • SSE solved by RK4 method = 29.8589927079. ‧. y. Nat. We have decrease the SSE from “55.3” to “29.8”. But it is still not good enough. Thus. io. sit. we will now add some perturbation to p1 and p2 to lower the SSE. After that, we add some. n. al. er. perturbation to the initial value as well. Repeating these two steps, we find an initial value and. i n U. v. a set of parameters which are relatively good. We can understand this part better by checking. Ch. engchi. out the Python code in appendix A.1. Result is as follow: • adjusted p1 = -0.00959911519495 • adjusted p2 = 0.232405855098 • initial value = (0.990787303889, 1.07654767318) • SSE solved by Heun’s method = 5.6053053617 • SSE solved by RK2 method = 5.58681233538 • SSE solved by RK4 method = 5.533732902. 28.

(39) The best result solved by this method is shown in Figure 4.7. The blue line is the original data of Kandelia’s Height, and the red line is our fitting curve.. 學. ‧ 國. 立. 政治大 Figure 4.7: Best result. ‧. In conclusion, we have decrease the final SSE to “5.533”. Compare to the result in former. n. Ch. y. 5.533 ≈ 0.16666 6.63916. sit. io. al. 1−. er. Nat. subsection SSE = 6.3916, we decrease the SSE by:. engchi. i n U. v. Which means we reduced about 16% of sum of square error. However we recall that in section 3.3, we build a model with SSE = “4.92591951”. It is still the best result until now. Thus we need to figure out some other method to decrease SSE.. 4.3.3 Using Analytical Solution Before moving on to the next section, we look at equation 4.3.2 again: Q ′ ( t ) = p1 Q2 ( t ) + p2 Q ( t ). We notice that this is the form of Bernoulli differential equation [2]. Thus, in this subsection, we will first use Matlab curve fitting tool to find suitable p1 and p2 . Then, we will find 29.

(40) analytical solution of this equation. Here is how it works:. Similar to subsection 4.3.1, we are going to fit the graph with Kandelia’s height as x axis and Kandelia’s growth velocity as y axis. Using the model below, we then find p1 and p2 by Matlab curve fitting tool. Q ′ ( t i ) = p1 Q2 ( t i ) + p2 Q ( t i ) Here, ( Q′ (ti ), Q(ti )), i = 0, · · · , n are the data of Kandelia’s velocity and Kandelia’s height. 政治大. at time ti . Suitable parameters we found are as follow. Also, we plot the result in Figure 4.8:. 立. • p1 = -0.008447. ‧. ‧ 國. 學. • p2 = 0.2117. n. er. io. sit. y. Nat. al. Ch. engchi. i n U. v. Figure 4.8: Graph of data and (4.3.1) We will now find the analytical solution of equation 4.3.2: Q ′ ( t ) = p1 Q2 ( t ) + p2 Q ( t ). ⇒ Q′ (t) + (− p2 ) Q(t) = p1 Q2 (t) ⇒ Q−2 (t) Q′ (t) + (− p2 ) Q−1 (t) = p1. 30.

(41) Here, we let w(t) = Q−1 (t), then we have w′ (t) = − Q−2 (t) Q′ (t). Substitute this into the above equation, we have the follow:. −w′ (t) + (− p2 )w(t) = p1 ⇒ w ′ ( t ) + p2 w ( t ) = − p1. (4.3.4). Equation 4.3.4 turned out to be a first order linear differential equation. Thus we solve it and substitute Q(t) back into the solution. We then have the following result:. Q(t) =. 立. p2 × 106 × e p2 t 6 e(c1 ×10 × p2 ) − p1 × 106 × e p2 t. 政治大. (4.3.5). Since p1 and p2 are already found. We now want to find suitable c1 in equation 4.3.5. Again,. ‧ 國. 學. using Matlab curve fitting tool, we have c1 = 5.717 × 10−5 . The result is shown in Figure 4.9. ‧. n. er. io. sit. y. Nat. al. Ch. engchi. i n U. v. Figure 4.9: Graph of data and (4.3.5) Overall, we have SSE = “5.646539” here. The result is almost the same as previous subsection SSE = “5.533733”. This shows that the solution found by genetic algorithm is similar to the analytical solution.. 31.

(42) 4.4. Model 4. Until now, the best result in section 4.3 is SSE = ‘5.533733”. However, we recall the result in section 3.3. There, SSE = “4.92591”, which is still better. Therefore, in this section, we give equation 4.3.2 a furthermore correction. We now give a closer examination on Figure 4.8. Starting from height in range 0 to 9, the data is higher than the fitting function. Next , in the range 10 to 20, the data is lower. After that, form range 21 to 25, it goes higher again. This difference seems like the behavior of a periodic factor. Thus, we build our new model by adding a sin term to equation 4.3.2. It is shown as. 政治大 Q (t) = p Q (t) + p Q(t) + asin(bQ(t) + c) 立. follow:. ′. 1. 2. 2. (4.4.1). ‧ 國. 學. Again, our goal is to find suitable p1 , p2 , a, b and c of equation 4.4.1. In subsection 4.4.1, we will use finite difference method. This substitutes term of Q′ (t) by Q(t) in our model. ‧. function. Then we can fit it by Table 2.3. Next, in subsection 4.4.2, we will fit equation 4.4.1. sit. y. Nat. by genetic algorithm directly. Notice that Table 4.3 is used here to fit our model. At last, in. io. al. n. similar to subsection 4.4.2.. er. subsection 4.4.3, we will use Matlab curve fitting tool to do the job. We will see the result is. Ch. e. ngch 4.4.1 Using Finite Difference Method. i. i n U. v. In this subsection, our goal is to find suitable parameters of equation 4.4.1: Q′ (t) = p1 Q2 (t) + p2 Q(t) + asin(bQ(t) + c). Here is how we do it. We first use the central difference method as our key ideal. That is, we will substitutes term of Q′ (t) by Q(t) in the above model equation. The equation will turn out to be a function with Q(t) and t. Thus, with some initial parameters, we can fit it by data of Kandelia’s height. Secondly, we use genetic algorithm to improve our model. That is, we give some perturbation to the initial value. Then, we compare the outcome depending on different initial value. 32.

(43) After that, we choose the best one. Finally, we use genetic algorithm again. This time, we give some perturbation on parameters. With these new parameters, we give initial value some perturbation. Repeating these two process, we will at last get a set of parameters depend on an initial value.We will now give a further detail of how this all works:. Let Q(ti ) be the data of Kandelia’s height at time ti , i = 1, · · · , n − 1. First, we substitute Q′ (t) in equation 4.4.1 by. Q(ti+1 )− Q(ti−1 ) t i +1 − t i −1. 政治大. Q ( t i +1 ) − Q ( t i −1 ) = p1 Q2 (ti ) + p2 Q(ti ) + asin(bQ(ti ) + c) t i +1 − t i −1. 立. ⇒ Q(ti+1 ) = Q(ti−1 ) + (ti+1 − ti−1 )( p1 Q2 (ti ) + p2 Q(ti ) + asin(bQ(ti ) + c). ‧ 國. 學. (4.4.2). ‧. We first give a set of initial coefficient by referring the result in subsection 4.3.2:. sit. y. Nat. • p1 = −0.00844749. n. al. er. io. • p2 = 0.21174495 • a = 0.0 • b = 0.0. Ch. engchi. i n U. v. • c = 0.0 Notice that these parameters can be chosen randomly. However, if the parameter is not chosen properly, the method by genetic algorithm may not work well. We will give a further discussion of this problem later. Anyway, by using parameters above we can solve equation 4.4.1 by Heun’s method, RK2 method and RK4 method. • SSE solved by Heun’s method = 55.7242229592 • SSE solved by RK2 method = 55.6657350223 • SSE solved by RK4 method = 55.3013536756 33.

(44) From our experience before, the result is always not so good at the first step. In order to lower SSE, we add some perturbation to the initial value. Thus we can get different data by using equation 4.4.2 within different initial value. We then use our new initial value to solve equation 4.4.2. The result is as follow : • adjusted initial value = (1.20182153439, 1.32129224963) • SSE solved by Heun’s method = 5.79278605472 • SSE solved by RK2 method = 5.7761112128. 政治大. • SSE solved by RK4 method = 5.7418495602. 立. The above result is pretty good. To lower SSE even more, we will now add some perturbation. ‧ 國. 學. to p1 , p2 , a, b and c. After that, we add some perturbation to the initial value as well. Repeating these two steps, we hope to find best parameters and the initial value they depends on. We have. ‧. the following result :. y. Nat. er. io. sit. • adjusted p1 = -0.00951510827755 • adjusted p2 = 0.230182300243. n. al. Ch. • adjusted a = -0.0262259316043. engchi. i n U. • adjusted b = -0.0147144948651 • adjusted c = -0.00726659819372 • initial value = (1.01627076858, 1.10278819877) • SSE solved by Heun’s method = 5.58823184026 • SSE solved by RK2 method = 5.56985208097 • SSE solved by RK4 method = 5.51999849112. 34. v.

(45) The best result solved by this method is shown in Figure 4.10. The blue line is our original data of Kandelia’s height, and the red line is our fitting curve. Further detail of how this get from Python is in appendix A.3. 立. 政治大. ‧ 國. 學 Figure 4.10: Best result. ‧ sit. y. Nat. Until now our best result in this section is SSE = “5.519998”. It slightly better than the result. io. n. al. er. in section 4.3. There, SSE = “5.533733”. Compare these two result: 1−. Ch. i n U. 5.519998 ≈ 0.002482 5.533733. engchi. v. This means, we only decrease about 0.24% of SSE. Also, recall the result in section 3.3,we have SSE = “4.9259”. It is still the best model until now. Therefore, we will try to improve it by using different method.. 4.4.2 Using Genetic Algorithm In this subsection, we again try to find parameters that fits the equation below: Q′ (t) = p1 Q2 (t) + p2 Q(t) + asin(bQ(t) + c) Different from the previous subsection, we will not change the term of Q′ (t) here. Instead, 35.

(46) we will use the data of Kandelia’s growth velocity as Q′ (t). Therefore, using Table 4.3, we can plot a graph with Kandelia’s height as x axis and Kandelia’s growth velocity as y axis. We then use the method of genetic algorithm to find suitable parameters. Here is how we do it:. First, we give a set of initial parameters. Notice that these parameters are chosen in purpose. We can see that in Figure 4.11 the equation generated by them won’t fit so badly. • p1 = −0.01 • p2 = 0.2 • a = −0.5. 立. • b = −0.3. ‧. ‧ 國. 學. • c = −0.1. 政治大. Here, the blue line represents the data of Kandelia’s growth rate, and the red line is our fitting. n. al. er. io. sit. y. Nat. curve.. Ch. engchi. i n U. v. Figure 4.11: Fit the data by genetic algorithm Next, in order to fit our data better, we add some perturbation to these parameters. Then, we compare the outcome of equation generated by different parameters. At last, we get the best result as follow. Also, we plot it out in Figure 4.12 36.

(47) • adjusted p1 = -0.0114876114475 • adjusted p2 = 0.262698012377 • adjusted a = -0.53344134757 • adjusted b = -0.318721070024 • adjusted c = -0.120528934405 Again, the blue line represents the data of Kandelia’s growth rate, and the red line is our fitting curve.. 立. 政治大. ‧. ‧ 國. 學. n. er. io. sit. y. Nat. al. Ch. engchi. i n U. v. Figure 4.12: Fit the data by genetic algorithm with parameters adjusted Finally, we solve equation 4.4.1 by Heun’s method, RK2 method and RK4 method. After that, we give initial value some perturbation as well. This can help us fit the original data better. A further detail of python code used in this subsection is in appendix A.4. • initial value for Heun’s method= 0.0489081076939 • initial value for RK2 method= 0.0489704328867 • initial value for RK4 method= 0.0486656124813 • SSE solved by Heun’s method = 1.68479162933 37.

(48) • SSE solved by RK2 method = 1.68477818706 • SSE solved by RK4 method = 1.68491878659 Final result is in Figure 4.13. Here, the blue line represents the data of Kandelia’s height, and the red line is our fitting curve.. 立. 政治大. ‧. ‧ 國. 學. Nat. er. io. sit. y. Figure 4.13: Best result. In conclusion, we have our best SSE = “1.68477818706”. This is the best result until now.. n. al. Ch. Thus, we compare this to SSE in section 3.3: 1−. engchi. i n U. v. 1.684778 ≈ 0.657975 4.9259. This means, we decrease about 65% of SSE. Although we are satisfied with the SSE we get, we still have to check this model does not over fit the data. Thus, we solve equation 4.4.1 with RK4 method. But this time, we doubled the range of time and double the grids as well. The result is as follow:. 38.

(49) 政治大. Figure 4.14: Predict of Kandelia’s height by equation 4.4.1. 立. ‧ 國. 學. In Figure 4.14, the blue line represents the original data of Kandelia’s height, and the red line is our fitting curve. While t > 30, the fitting curve grows slowly and approaches stable.. er. io. sit. Nat. 4.4.3 Using Matlab Curve Fitting Tool. y. ‧. This is more reasonable to the real world. Therefore, we take this model as our final result.. In this section, we will compare the result from genetic algorithm to the result using Matlab. al. n. v i n curve fitting tool. We want to show Thus, we again set equation 4.4.1 as C hthat they are similar. engchi U our model: Q′ (t) = p1 Q2 (t) + p2 Q(t) + asin(bQ(t) + c). Again, we use the data in Table 4.3 to construct a graph. Our goal is to find suitable parameters by Matlab curve fitting tool. The corresponding parameters are as follow. Also, we have our result shown in Figure 4.15: • p1 = -0.011 • p2 = 0.2601 • a = -0.5341 • b = -0.3017 39.

(50) • c = -0.1563. 治政大 Figure 4.15: Graph of data and equation(4.4.1) 立 ‧ 國. 學. Using the parameters above, we can find a numerical solution of equation4.4.1 via Heun’s method, RK2 method and RK4 method. Also we give the initial value some perturbation in. ‧. order to decrease the SSE. Results are as below:. y. Nat. n. al. er. io. • initial value for RK2 method= -0.0153907948845. sit. • initial value for Heun’s method= -0.0151605021498. i n C • initial value for RK4 method=h -0.0152950386194 engchi U. v. • SSE solved by Heun’s method = 2.12606450927 • SSE solved by RK2 method = 2.12599915696 • SSE solved by RK4 method = 2.12601296975 The best result solved by this method is shown in Figure 4.16. Here, the blue line represents the data of Kandelia’s height, and the red line is our fitting curve.. 40.

(51) 立. 政治大. Figure 4.16: Best result. ‧ 國. 學. Overall, we have SSE = “2.215999” . It is slightly worse than the result in previous subsection (SSE = “1.68477818706”). However, the result here is the second best in this paper. Also,. n. al. er. io. sit. y. Nat. subsection 4.4.2.. ‧. we notice that parameters we found by Matlab curve fitting tool is similar to those we found in. Ch. engchi. 41. i n U. v.

(52) Chapter 5 Conclusion 政治大 We sort out all the mathematical models in table 5.1. Here is a quick review of what we have 立. ‧ 國. 學. done in this paper:. In chapter 3, we built two model functions. In section 3.2 the model is constructed by logistic. ‧. curve. There, SSE = “5.38952958”. Then, in section 3.3 we improve this model by adding a sin function. By doing so, we reduced SSE to “4.92591951”. It means, we decrease about 8.6% of. io. sit. y. Nat. first SSE.. n. al. er. Next in chapter 4, we built two model of differential equation. In section 4.3, we construct. i n U. v. a model from differential form of logistic curve. There, we have outcome SSE = “5.533733”,. Ch. engchi. which is not a good one. However, in section 4.4, we improve the model equation by adding a sin function. Then, we reduced SSE to “1.68477819”. That is, we decrease about 65% of SSE in section 3.3. In conclusion, we find out the result is the best with following three condition: • fit data: ( Q′ (t), Q(t)) • fit model: Q′ (t) = p1 Q2 (t) + p2 Q(t) + asin(bQ(ti ) + c) • genetic algorithm as our fit method However, this method also has its disadvantages. Like the Matlab curve fitting tool, if the model function with initial parameters has a large SSE with original data. Then the parameters. 42.

(53) may probably not converge. Therefore, how to find a suitable initial parameters may be a good question for further research.. Table 5.1: Conclusion Section Section 3.2 Section 3.3. Fit model Q(t) = 1+ AeB−Bkt Q(t) = 1+ AeB−Bkt + asin(bt + c). Section 4.3.1. Q ′ ( t ) = p1 Q2 ( t ) + p2 Q ( t ). 政治大. 立. ‧. ‧ 國. = p1 Q2 ( t i ) + p2 Q ( t i ). Q ′ ( t ) = p1 Q2 ( t ) + p2 Q ( t ). n. al. er. io. sit. y. Nat. Section 4.3.3. Q(ti+1 )− Q(ti−1 ) t i +1 − t i −1. 學. Section 4.3.2. Fit Method Matalb curve fitting tool Matalb curve fitting tool Heun’s Regression RK2 RK4 Heun’s Finite difference RK2 RK4 Heun’s Initial value adjust RK2 RK4 Heun’s Coefficient adjust RK2 RK4 Analytical Solution Heun’s Finite difference RK2 RK4 Heun’s Initial value adjust RK2 RK4 Heun’s Coefficient adjust RK2 RK4 Heun’s Genetic algorithm RK2 RK4 Heun’s Curve fitting tool RK2 RK4. i n U. Section 4.4.1. Q(ti+1 )− Q(ti−1 ) = t i +1 − t i −1 2 p1 Q ( t i ) + p2 Q ( t i ) +. Section 4.4.2. Q′ (t) = p1 Q2 (t) + p2 Q(t) + asin(bQ(ti ) + c). Section 4.4.3. Q′ (t) = p1 Q2 (t) + p2 Q(t) + asin(bQ(t) + c). Ch. e nasin g(cbQh(tii ) + c). 43. v. SSE 5.38952958 4.92591951 6.39167107 6.39166484 6.39165420 55.7242174 55.6657295 55.3013481 30.2074083 30.1259883 29.8589927 5.605305 5.586812 5.533733 5.646539 55.724230 55.665735 55.301354 5.7927861 5.7761112 5.7418496 5.58823184 5.56985208 5.51999849 1.68479163 1.68477819 1.68491879 2.12606451 2.12599916 2.12601297.

(54) Appendix A Code use in paper. ‧. %pylab inline. y. sit er. al. n. In[2]:. io. Heun’s Method. Nat. Some Tools. 1. 學. In[1]: 1. 立. Code 01. ‧ 國. A.1. 政治大. Ch. engchi. def ode_Heun(f,t0,tf,y0,y1,n):. 2. t = np.linspace(t0,tf,n). 3. y = list([y0,y1]). 4. for i in range(1,n-1):. i n U. 5. h = t[i+1]-t[i]. 6. fk = f(t[i],y[-1]). 7. fk1 = f(t[i+1],y[-1]+h*fk). 8. y.append(y[-1]+h*(fk+fk1)/2.0). 9. return y RK2 Method. 44. v.

(55) In[3]: 1. def ode_RK2(f,t0,tf,y0,y1,n):. 2. t = linspace(t0,tf,n). 3. y = list([y0,y1]). 4. for i in range(1,n-1):. 5. h = t[i+1]-t[i]. 6. k1 = h*f(t[i],y[-1]). 7. k2 = h*f(t[i]+h/2.0,y[-1]+k1/2.0). 8. y.append(y[-1]+k2) return y. 9. 政治大. 立. RK4 Method. def ode_RK4(f,t0,tf,y0,y1,n):. 2. t = linspace(t0,tf,n). 3. y = list([y0,y1]). 4. for i in range(1,n-1):. sit. n. er. io. al. y. Nat. 1. ‧. ‧ 國. 學. In[4]:. i n U. v. 5. h = t[i+1]-t[i]. 6. k1 = h*f(t[i],y[-1]). 7. k2 = h*f(t[i]+h/2.0,y[-1]+k1/2.0). 8. k3 = h*f(t[i]+h/2.0,y[-1]+k2/2.0). 9. k4 = h*f(t[i]+h,y[-1]+k3). Ch. engchi. y.append(y[-1]+(k1+2*k2+2*k3+k4)/6). 10. return y. 11. SSE In[5]: 1. def SSE_finitedifference(y_estimate,y_real):. 2. SSE_Temp = 0. 3. for i in range(len(y_real)): 45.

(56) SSE_Temp = SSE_Temp + (y_real[i]-y_estimate[i])**2. 4. 5. SSE_Temp = sqrt(SSE_Temp). 6. return SSE_Temp Data In[6]:. 1. y =. [0.0,0.1,0.3,0.6,0.8,1,2.5,3,4.5,5.7,7.9,9.2,11.5,12,13.2,14.2,15,5.7,1. 政治大 [0.15,0.25,0.25,0.2,0.85,1.0,1.0,1.35,1.7,1.75,1.8,1.4,0.85,1.1,0.9,0.7 立. 3. n = len(y). 4. n_cdy = len(cdy). 5. t = range(n). ‧ er. io. sit. y. Nat. In[7]:. ‧ 國. cdy =. 學. 2. 1. xlabel(’t:Time (weeks)’). 2. ylabel(’u(t):Height (cm)’). 3. plot(t,y,marker = ’o’). n. al. Ch. engchi. i n U. v. 1. Finite Difference Method 求 y′ = αy2 + βy 的最佳 α 和 β 利⽤ central difference method 46.

(57) y k +1 − y k −1 = αy2k + βyk t k +1 − t k −1. → yk+2 = yk + 2αy2k+1 + 2βyk+1 Step1. 利⽤最⼩平⽅法求出最佳 α 和 β In[8]: 1. A = zeros((n-2,2)). 2. b = zeros(n-2). 3. 立. for i in range(n-2):. 學. 5. A[i,0] = y[i+1]**2. 6. A[i,1] = y[i+1]. 7. b[i] = 0.5*(y[i+2] - y[i]). 11. 12. y. sit. b = mat(b). al. er. 10. n. A = mat(A). io. 9. Nat. 8. ‧. ‧ 國. 4. 政治大. Ch. engchi. i n U. v. z = linalg.solve(A.T*A,A.T*b.T). 13. 14. alpha = float(z[0]). 15. beta = float(z[1]). 16. z matrix([[-0.00844749],[ 0.21174495]]). Step.2 接下來利⽤上⾯求出的 α 和 β，搭配上數值偏微分⽅法來得到⼀組數據逼近原始資料也就是說 Given {(yi , ti )}in=1 ，我們要⽤ Heun’s Method、RK2 Method 和 RK4 Method 47.

(58) 去找到⼀組資料來逼近原始資料 Settintg up Data 如上⾯給定的 In[9]: 1. t0 = t[0]. 2. tf = t[-1]. 3. y0 = 0.0. 4. y1 = 0.1. 5. 6. #要解的微分⽅程. 立. 7. ‧ 國. 9. 學. 8. 政治大. def cdy_fcn(t,y): return alpha*y**2+beta*y. ‧. Heun’s Method. io. sit. y. Nat. In[10]:. y_Heun = ode_Heun(lambda x,y: cdy_fcn(x,y),t0,tf,y0,y1,n). 2. SSE_Heun0 = SSE_finitedifference(y_Heun,y). 3. print SSE_Heun0. n. al. er. 1. Ch. engchi. i n U. v. 55.7242174456. RK2 Method In[11]: 1. y_RK2 = ode_RK2(lambda x,y: cdy_fcn(x,y),t0,tf,y0,y1,n). 2. SSE_RK2a = SSE_finitedifference(y_RK2,y). 3. print SSE_RK2a 55.66572951. 48.

(59) RK4 Method In[12]: 1. y_RK4 = ode_RK4(lambda x,y: cdy_fcn(x,y),t0,tf,y0,y1,n). 2. SSE_RK4a = SSE_finitedifference(y_RK4,y). 3. print SSE_RK4a 55.3013480532. Step3.** 對 y0 ，y1 做擾動，希望能降低 SSE. y k +2. 立. ‧ 國. y_1 = [0.0 , 0.1]. 2. for i in range(n-2):. Nat. al. er. io. sit. y. y_1.append(y_1[-2]+2*alpha*y_1[-1]**2+2*beta*y_1[-1]). 4. #最⼩平⽅差. 6. SSE_0 = SSE_finitedifference(y,y_1). 7. print SSE_0. n. 5. Ch. engchi. In[14]: 1. count = 0. 2. SSE = SSE_0. 3. while count <1000:. 4. k +1. ‧. 1. 3. 2 k +1. k. 學. In[13]:. 政治大 = y + 2αy + 2βy. y_temp = [0.0, 0.0]. 5. 6. #對起始點增加擾動. 7. e1 = random.uniform(-0.1,0.1). 8. e2 = random.uniform(-0.1,0.1) 49. i n U. v.

(60) 9. y_temp[0] = y_1[0]+e1. 10. y_temp[1] = y_1[0]+e2. 11. 12. # ⽣成新的Data. 13. for i in range(n-2): y_temp.append(y_temp[-2]+2*alpha*y_temp[-1]**2+2*beta*y_temp. 14. [-1]) 15. 16. #計算新的最⼩平⽅差. 17. SSE_temp = SSE_finitedifference(y,y_temp). 立. 學. #⽐較最⼩平⽅差. 20. if SSE_temp < SSE:. 21. SSE = SSE_temp. 22. y_1 = y_temp. sit. al. n. count += 1. io. 24. er. 23. y. Nat. 19. ‧. ‧ 國. 18. 政治大. 25. 26. print SSE. 27. print alpha,beta. 28. print y_1[0], y_1[1]. Ch. engchi. i n U. v. 5.659286856 -0.0084474869693 0.211744951986 1.16903645355 1.32198681127. Heun’s Method In[15]: 1. y_Heun = ode_Heun(lambda x,y: cdy_fcn(x,y),t0,tf,y_1[0],y_1[1],n) 50.

(61) 2. SSE_Heun0 = SSE_finitedifference(y_Heun,y). 3. print SSE_Heun0 30.2074082818. RK2 Method In[16]: 1. y_RK2 = ode_RK2(lambda x,y: cdy_fcn(x,y),t0,tf,y_1[0],y_1[1],n). 2. SSE_RK2a = SSE_finitedifference(y_RK2,y). 3. print SSE_RK2a. 立. 30.1259883441. ‧ 國. ‧. In[17]:. 學. RK4 Method. 政治大. Nat. 2. SSE_RK4a = SSE_finitedifference(y_RK4,y). 3. print SSE_RK4a. er. n. al. 29.8589927079. Ch. engchi. Step.4 對 α、β 和 y0 做擾動，降低 SSE 有了 Step.3 修正過的 y0 我們接著來修正 α 和 β In[18]: 1. flag = 0. 2. while flag < 100:. 3. count = 0. 4. SSE = SSE_0. 5. while count <1000:. 6. sit. y. y_RK4 = ode_RK4(lambda x,y: cdy_fcn(x,y),t0,tf,y_1[0],y_1[1],n). io. 1. y_temp = [y_1[0], y_1[1]] 51. i n U. v.

(62) 7. 8. #對起始點增加擾動. 9. e1 = random.uniform(-0.001,0.001). 10. e2 = random.uniform(-0.001,0.001). 11. alpha_temp = alpha + e1. 12. beta_temp = beta + e2. 13. 14. # ⽣成新的Data. 15. for i in range(n-2):. 16. 政治大. y_temp.append(y_temp[-2]+2*alpha_temp*y_temp[-1]**2+2*. 立. beta_temp*y_temp[-1]). ‧ 國. 學. 17. #計算新的最⼩平⽅差. 18. SSE_temp = SSE_finitedifference(y,y_temp). ‧. #⽐較最⼩平⽅差. 21. if SSE_temp < SSE:. 22. SSE = SSE_temp. 23. y_1 = y_temp. 24. alpha = alpha_temp. 25. beta = beta_temp. sit. io. n. al. y. Nat. 20. er. 19. Ch. engchi. 26. 27. count += 1. 28. 29. count = 0. 30. while count <1000:. 31. y_temp = [0.0, 0.0]. 32. 33. #對起始點增加擾動. 52. i n U. v.

(63) 34. e1 = random.uniform(-0.1,0.1). 35. e2 = random.uniform(-0.1,0.1). 36. y_temp[0] = y_1[0]+e1. 37. y_temp[1] = y_1[0]+e2. 38. 39. # ⽣成新的Data. 40. for i in range(n-2): y_temp.append(y_temp[-2]+2*alpha*y_temp[-1]**2+2*beta*. 41. y_temp[-1]) 42. 政治大. 立. #計算新的最⼩平⽅差. 44. SSE_temp = SSE_finitedifference(y,y_temp). 45. ‧. ‧ 國. 學. 43. #⽐較最⼩平⽅差. 47. if SSE_temp < SSE:. 48. SSE = SSE_temp. 49. y_1 = y_temp. sit. n. 50. Ch. count += 1. 51. 52. er. io. al. y. Nat. 46. engchi. flag+=1. 53. print SSE. 54. print alpha,beta. 55. print y_1[0], y_1[1] 5.39888930477 -0.00959911519495 0.232405855098 0.990787303889 1.07654767318. In[19]:. 53. i n U. v.

(64) 1. t0 = t[0]. 2. tf = t[-1]. 3. y0 = y_1[0]. 4. y1 = y_1[1]. 5. #要解的微分⽅程. 6. def cdy_fcn(t,y):. 7. return alpha*y**2+beta*y Heun’s Method. 政治大 y_Heun = ode_Heun(lambda x,y: cdy_fcn(x,y),t0,tf,y0,y1,n) 立 n[20]:. 1. SSE_Heun0 = SSE_finitedifference(y_Heun,y). 3. print SSE_Heun0. ‧ 國. y. sit er. io. al. v i n ode_RK2(lambda C x,y: h e ncdy_fcn(x,y),t0,tf,y0,y1,n) gchi U n. In[21]:. Nat. RK2’s Method. ‧. 5.6053053617. 學. 2. 1. y_RK2 =. 2. SSE_RK2a = SSE_finitedifference(y_RK2,y). 3. print SSE_RK2a 5.58681233538. RK4 Method In[22]: 1. y_RK4 = ode_RK4(lambda x,y: cdy_fcn(x,y),t0,tf,y0,y1,n). 2. SSE_RK4a = SSE_finitedifference(y_RK4,y). 3. print SSE_RK4a. 54.

(65) 5.533732902. A.2. Code 02. In[1]: 1. %pylab inline Some Tools Heun’s Method In[2]:. 立. def ode_Heun(f,t0,tf,y0,n):. 學. ‧ 國. t = np.linspace(t0,tf,n). 3. y = list([y0]). 4. for i in range(n-1):. ‧. 2. y. Nat. 6. fk = f(t[i],y[-1]). 7. fk1 = f(t[i+1],y[-1]+h*fk). 8. y.append(y[-1]+h*(fk+fk1)/2.0). n. al. sit. h = t[i+1]-t[i]. io. 5. Ch. engchi. 9. 10. return y RK2 Method In[3]:. 1. def ode_RK2(f,t0,tf,y0,n):. 2. t = linspace(t0,tf,n). 3. y = list([y0]). 4. for i in range(n-1):. 5. h = t[i+1]-t[i]. 6. k1 = h*f(t[i],y[-1]) 55. er. 1. 政治大. i n U. v.

(66) 7. k2 = h*f(t[i]+h/2.0,y[-1]+k1/2.0). 8. y.append(y[-1]+k2). 9. return y RK4 Method In[4]:. 1. def ode_RK4(f,t0,tf,y0,n):. 2. t = linspace(t0,tf,n). 3. y = list([y0]). 4. for i in range(n-1):. 立. 政治大. h = t[i+1]-t[i]. 6. k1 = h*f(t[i],y[-1]). 7. k2 = h*f(t[i]+h/2.0,y[-1]+k1/2.0). 8. k3 = h*f(t[i]+h/2.0,y[-1]+k2/2.0). 9. k4 = h*f(t[i]+h,y[-1]+k3). ‧ 國. y. sit. n. al. er. io. return y SSE In[5]:. 1. Ch. engchi. i n U. v. def SSE_cdy(y_estimate,y_real):. 2. SSE_Temp = 0. 3. for j in range(len(y_real)): SSE_Temp = SSE_Temp + (y_estimate[j]-y_real[j])**2. 4. 5. SSE_Temp = sqrt(SSE_Temp). 6. return SSE_Temp In[6]:. 1. 2. ‧. Nat. y.append(y[-1]+(k1+2*k2+2*k3+k4)/6). 10. 11. 學. 5. def SSE_numerical(y_estimate,y_real,gap): y_T = [] 56.

(67) 3. SSE_Temp = 0.0. 4. for i in range(len(y_real)):. 5. y_T.append(y_estimate[i*gap]). 6. SSE_Temp = SSE_Temp+(y_T[i]-y_real[i])**2. 7. SSE_Temp = sqrt(SSE_Temp). 8. return SSE_Temp Data In[7]:. 政治大 [0.0,0.1,0.3,0.6,0.8,1,2.5,3,4.5,5.7,7.9,9.2,11.5,12,13.2,14.2,15,15.7, 立. 2. cdy =. ‧ 國. y =. 學. 1. ‧. [0.15,0.25,0.25,0.2,0.85,1.0,1.0,1.35,1.7,1.75,1.8,1.4,0.85,1.1,0.9,0.7. 5. t = range(n). sit. n_cdy = len(cdy). n. al. er. 4. io. n = len(y). y. Nat. 3. In[8]:. Ch. engchi. 1. xlabel(’t:Time (weeks)’). 2. ylabel(’u(t):Height (cm)’). 3. plot(t,y,marker = ’o’). 57. i n U. v.

(68) 2. 利⽤ (y, cdy) 來找出 y′ = αy2 + βy 的最佳 α 和 β. 政治大. Step1. 利⽤最⼩平⽅法求出最佳 α 和 β In[9]:. 學. y_5 =. ‧ 國. 1. 立. [0.0,0.1,0.3,0.6,0.8,1,2.5,3,4.5,5.7,7.9,9.2,11.5,12,13.2,14.2,15,15.7,. ‧. xlabel(’u(t):Height (cm)’). 3. ylabel(’u(t):velocity (cm/(weeks))’). 4. plot(y_5,cdy,marker = ’o’,color = ’b’,). n. Ch. engchi. In[10]: 1. # 利⽤最⼩平⽅法求出 alpha 和 beta. 2. A = zeros((n_cdy,2)). 3. b = zeros(n_cdy) 58. sit er. io. al. y. Nat. 2. i n U. v.

(69) 4. 5. for i in range(n_cdy):. 6. A[i,0] = y[i]**2. 7. A[i,1] = y[i]. 8. b[i] = cdy[i]. 9. 10. A = mat(A). 11. b = mat(b). 12. 13. 政治大. z = linalg.solve(A.T*A,A.T*b.T). 立. 學. ‧ 國. 14. alpha = float(z[0]). 16. beta = float(z[1]). 17. z. ‧. 15. y. Nat. n. al. 1. 2. 3. er. io. In[11]:. sit. matrix([[-0.00863605],[ 0.21215756]]). Ch. engchi. i n U. v. def cdy_fcn(t,y): return alpha*y**2+beta*y. 4. 5. cdy_estimate = []. 6. for i in range(n_cdy):. 7. cdy_estimate.append(cdy_fcn(t,y_5[i])). 8. 9. xlabel(’u(t):Height (cm)’). 10. ylabel(’u(t):velocity (cm/(weeks))’). 11. plot(y_5,cdy,marker = ’o’,color = ’b’,) 59.