1 00:00:04,550 --> 00:00:06,900 CLIVE MOLER: The Lotka-Volterra Altera predator 2 00:00:06,900 --> 00:00:12,820 prey equations are the granddaddy of all models 3 00:00:12,820 --> 00:00:16,950 involvement competition between species. 4 00:00:16,950 --> 00:00:19,480 They are the foundation of fields 5 00:00:19,480 --> 00:00:22,860 like mathematical ecology. 6 00:00:22,860 --> 00:00:29,440 Think of the two species as rabbits and foxes 7 00:00:29,440 --> 00:00:36,850 or moose and wolves or little fish in big fish. 8 00:00:36,850 --> 00:00:42,290 Y1 represents the prey, who would live peacefully 9 00:00:42,290 --> 00:00:46,020 by themselves if there were no predators. 10 00:00:46,020 --> 00:00:49,920 I've chosen the units of time and population 11 00:00:49,920 --> 00:00:51,840 so that the coefficients in front 12 00:00:51,840 --> 00:00:57,130 of the leading linear terms are one. 13 00:00:57,130 --> 00:01:04,370 So y1 prime equals y1 represents exponential growth 14 00:01:04,370 --> 00:01:07,625 of the prey in the absence of any predators. 15 00:01:10,610 --> 00:01:15,460 The predators need the prey, live on the prey. 16 00:01:15,460 --> 00:01:19,500 So in the absence of any prey, this minus sign 17 00:01:19,500 --> 00:01:21,090 is all important. 18 00:01:21,090 --> 00:01:27,930 So y2 prime equals minus y2 represents exponential decay. 19 00:01:27,930 --> 00:01:31,900 And the predators die off exponentially 20 00:01:31,900 --> 00:01:36,970 in the absence of any prey. 21 00:01:36,970 --> 00:01:39,680 But then here are the non-linear terms. 22 00:01:39,680 --> 00:01:45,300 These are like logistics terms, except with the interaction 23 00:01:45,300 --> 00:01:47,580 between the two species. 24 00:01:47,580 --> 00:01:52,880 The growth of Y1 is limited by the presence of y2. 25 00:01:52,880 --> 00:02:04,390 So y1 will grow until this term becomes 0 one y2 read reaches 26 00:02:04,390 --> 00:02:06,440 mu2. 27 00:02:06,440 --> 00:02:13,510 On the other hand, the decay of y2 28 00:02:13,510 --> 00:02:18,000 becomes 0 when y1 reaches mu1. 29 00:02:21,430 --> 00:02:25,940 To complete this specification, we need the initial conditions. 30 00:02:25,940 --> 00:02:30,220 So we have two values eta 1 and eta 2, 31 00:02:30,220 --> 00:02:35,150 which are the initial values of y1 and y2. 32 00:02:35,150 --> 00:02:40,760 So these four parameters, two mus and two etas, 33 00:02:40,760 --> 00:02:45,885 are the four parameters we have in our predator prey model. 34 00:02:49,340 --> 00:02:52,820 Don't worry about the fact that these are continuous variables 35 00:02:52,820 --> 00:02:57,110 and that we can have non-integer numbers of individuals. 36 00:02:57,110 --> 00:03:02,580 We can have half of a rabbit or a tenth of a moose. 37 00:03:02,580 --> 00:03:05,200 These are, after all, models that 38 00:03:05,200 --> 00:03:09,280 are idealized versions of what's happening in nature. 39 00:03:14,160 --> 00:03:19,740 The critical points are when the derivatives become 0. 40 00:03:19,740 --> 00:03:21,920 There's a critical point at the origin. 41 00:03:21,920 --> 00:03:26,140 But the interesting one is when these terms become 0. 42 00:03:29,920 --> 00:03:32,800 So that's the point where y1 and y2 43 00:03:32,800 --> 00:03:35,350 are equal to the mu1 and mu2. 44 00:03:38,510 --> 00:03:41,380 We have to look at the Jacobian. 45 00:03:41,380 --> 00:03:44,420 And here's the Jacobian in general. 46 00:03:44,420 --> 00:03:49,650 And at the critical point, the Jacobian-- here's the Jacobian. 47 00:03:49,650 --> 00:03:55,810 The eigenvalues of the Jacobian are plus or minus I. 48 00:03:55,810 --> 00:04:02,730 And so this critical point is a stable center 49 00:04:02,730 --> 00:04:05,035 with a period 2pi. 50 00:04:11,440 --> 00:04:13,790 So these are non-linear equations. 51 00:04:13,790 --> 00:04:16,360 We can't express the solution in terms 52 00:04:16,360 --> 00:04:19,640 of simple analytic functions. 53 00:04:19,640 --> 00:04:22,990 We have to compute them numerically. 54 00:04:22,990 --> 00:04:26,370 But we do know this about their behavior. 55 00:04:26,370 --> 00:04:30,820 If the initial conditions are close to the critical point, 56 00:04:30,820 --> 00:04:33,410 the solution is periodic. 57 00:04:33,410 --> 00:04:36,740 The period is close to 2pi. 58 00:04:36,740 --> 00:04:41,310 And the orbit is close to an ellipse. 59 00:04:41,310 --> 00:04:45,220 On the other hand, if the initial conditions are far 60 00:04:45,220 --> 00:04:48,060 from the critical point, then it turns out 61 00:04:48,060 --> 00:04:51,210 the solution is still periodic. 62 00:04:51,210 --> 00:04:56,660 But the period is greater than 2pi and the orbit 63 00:04:56,660 --> 00:05:00,060 is far from an ellipse. 64 00:05:04,680 --> 00:05:08,990 OK, let's bring up MATLAB and compute a solution. 65 00:05:08,990 --> 00:05:10,450 We need the parameters. 66 00:05:10,450 --> 00:05:15,410 Here's a mu and an eta. 67 00:05:15,410 --> 00:05:18,830 And I'll set the differential equation. 68 00:05:18,830 --> 00:05:26,410 And then choose ODE 45 and we'll integrate from 0 to 25. 69 00:05:26,410 --> 00:05:28,340 And here's the solution. 70 00:05:28,340 --> 00:05:30,080 It's periodic. 71 00:05:30,080 --> 00:05:32,390 A predator and a prey. 72 00:05:32,390 --> 00:05:36,550 And it looks like the period, it's 73 00:05:36,550 --> 00:05:41,120 returning back to the initial condition is 100 and 400. 74 00:05:41,120 --> 00:05:45,420 And it's returning back to that along about here. 75 00:05:45,420 --> 00:05:47,420 And we've integrated over something 76 00:05:47,420 --> 00:05:49,960 more than three periods. 77 00:05:49,960 --> 00:05:56,900 I happen to know that the period is 6.5. 78 00:05:56,900 --> 00:06:02,860 And so I can-- I'm going to set-- I want 79 00:06:02,860 --> 00:06:05,570 to compute to higher accuracy. 80 00:06:05,570 --> 00:06:08,940 10 to the minus 6. 81 00:06:08,940 --> 00:06:13,650 And let's capture the solution and integrate 82 00:06:13,650 --> 00:06:14,840 over three periods. 83 00:06:23,810 --> 00:06:29,460 It generated 337 points. 84 00:06:29,460 --> 00:06:37,730 Let's plot it with the finer dots. 85 00:06:37,730 --> 00:06:42,610 And we can see here, we've gone over three periods 86 00:06:42,610 --> 00:06:49,310 and return back to our initial values of 100 and 300. 87 00:06:49,310 --> 00:06:56,260 And now I'm going to use something 88 00:06:56,260 --> 00:07:00,290 that'll show off the periodicity of function 89 00:07:00,290 --> 00:07:01,940 in MATLAB called Comet. 90 00:07:15,350 --> 00:07:15,850 One orbit. 91 00:07:18,760 --> 00:07:21,520 Two orbits. 92 00:07:21,520 --> 00:07:22,265 Three orbits. 93 00:07:29,210 --> 00:07:32,150 Determining the period of a periodic solution 94 00:07:32,150 --> 00:07:36,530 is often the important part of a calculation. 95 00:07:36,530 --> 00:07:42,190 In the MATLAB ODE suite, this is done with an event handler. 96 00:07:42,190 --> 00:07:51,790 So I'm going to use ODE's set to provide an event handler called 97 00:07:51,790 --> 00:07:53,260 Pit Stop. 98 00:07:53,260 --> 00:07:57,960 And here's the code in Pit Stop. 99 00:07:57,960 --> 00:08:06,630 This code is called every time step in the integration. 100 00:08:06,630 --> 00:08:10,310 And I'm not going to go into detail here, 101 00:08:10,310 --> 00:08:15,160 but it computes a function g that we're 102 00:08:15,160 --> 00:08:17,530 looking to see when this is 0. 103 00:08:17,530 --> 00:08:22,660 And when g returns to 0, the integration is stopped. 104 00:08:22,660 --> 00:08:28,380 I'm going to-- I have a display function in here. 105 00:08:28,380 --> 00:08:29,990 Ordinarily, this wouldn't be here. 106 00:08:29,990 --> 00:08:31,760 But I want to look at this to see 107 00:08:31,760 --> 00:08:39,030 how the calculation is done. 108 00:08:39,030 --> 00:08:41,760 It says, terminate when y returns to the point 109 00:08:41,760 --> 00:08:45,770 where its angle mu is the same as the angle between eta 110 00:08:45,770 --> 00:08:46,990 and mu. 111 00:08:46,990 --> 00:08:52,440 This is more reliable than just finding 112 00:08:52,440 --> 00:08:59,500 when the solution returns back to its initial condition. 113 00:08:59,500 --> 00:09:03,820 So let's run ODE 45 and tell it to integrate over 114 00:09:03,820 --> 00:09:06,922 an infinite time step, over an infinite interval. 115 00:09:06,922 --> 00:09:08,630 That's not going to happen, because we're 116 00:09:08,630 --> 00:09:13,710 going to stop as soon as g is 0 but with the option vector set 117 00:09:13,710 --> 00:09:16,260 with this event handler. 118 00:09:16,260 --> 00:09:17,830 So here it is. 119 00:09:17,830 --> 00:09:23,350 And here's the output produced by Pit Stop 120 00:09:23,350 --> 00:09:25,630 as we integrate along. 121 00:09:25,630 --> 00:09:30,130 Here's the values that it's found. 122 00:09:30,130 --> 00:09:36,090 And here's where g returns to 0. 123 00:09:36,090 --> 00:09:41,940 I've produced a t vector with 117 values in it. 124 00:09:41,940 --> 00:09:49,310 And the final value of t is this value 125 00:09:49,310 --> 00:09:54,530 that I said was the period of this solution. 126 00:09:54,530 --> 00:09:57,720 So that's how the period was determined 127 00:09:57,720 --> 00:10:02,805 by this event handling feature in the ODE solvers. 128 00:10:06,110 --> 00:10:08,410 Great. 129 00:10:08,410 --> 00:10:12,260 I have a program called Predator Prey that's 130 00:10:12,260 --> 00:10:18,820 in the collection of programs that comes with NCM, Numerical 131 00:10:18,820 --> 00:10:21,160 Computing with MATLAB. 132 00:10:21,160 --> 00:10:25,375 My book that's available on the MathWorks website. 133 00:10:29,610 --> 00:10:38,830 Predator prey offers this graphic user interface 134 00:10:38,830 --> 00:10:41,810 to demonstrate what we've been talking about the predator prey 135 00:10:41,810 --> 00:10:44,660 equations. 136 00:10:44,660 --> 00:10:48,950 The top display shows the phase plane plot, the plot 137 00:10:48,950 --> 00:10:52,690 of prey versus predator. 138 00:10:52,690 --> 00:10:58,840 And the bottom display shows the time series plot, the plot 139 00:10:58,840 --> 00:11:00,870 of the two populations. 140 00:11:00,870 --> 00:11:02,920 And initially, it's set at the conditions 141 00:11:02,920 --> 00:11:04,880 we've been talking about. 142 00:11:04,880 --> 00:11:12,540 Here's 400 rabbits and 100 foxes around the critical point 143 00:11:12,540 --> 00:11:17,370 of 300 rabbits and 200 foxes. 144 00:11:17,370 --> 00:11:22,110 And here's the period of 6.5 some odd 145 00:11:22,110 --> 00:11:24,250 that we've been working with. 146 00:11:24,250 --> 00:11:27,510 Now, it says drag either dot. 147 00:11:27,510 --> 00:11:34,530 And you can change the initial conditions 148 00:11:34,530 --> 00:11:37,380 or the critical point. 149 00:11:37,380 --> 00:11:40,280 If I bring the initial conditions close 150 00:11:40,280 --> 00:11:46,200 to the critical point, then the phase plane plot 151 00:11:46,200 --> 00:11:48,480 becomes close to ellipse. 152 00:11:48,480 --> 00:11:53,760 And the period becomes close to 2pi. 153 00:11:53,760 --> 00:12:00,590 This is 6.28, which is twice 3.14. 154 00:12:00,590 --> 00:12:08,390 But if I change it so that the two are far away 155 00:12:08,390 --> 00:12:13,310 from each other, then the phase plane plot 156 00:12:13,310 --> 00:12:21,250 becomes quite different from an ellipse. 157 00:12:21,250 --> 00:12:22,550 It's always periodic. 158 00:12:22,550 --> 00:12:28,400 That's an amazing thing about these nonlinear equations. 159 00:12:28,400 --> 00:12:32,090 They always have a periodic solution. 160 00:12:32,090 --> 00:12:38,490 But now as you can see, the period is greater than 2pi. 161 00:12:38,490 --> 00:12:47,340 And the solutions are very far from sinusoids. 162 00:12:47,340 --> 00:12:51,860 So that's the pred prey app that's 163 00:12:51,860 --> 00:12:59,520 available with the programs in from the mathworks website 164 00:12:59,520 --> 00:13:04,705 under NCM for Numerical Computing with MATLAB. 165 00:13:10,910 --> 00:13:13,780 Whoops, I stand corrected. 166 00:13:13,780 --> 00:13:20,250 If you Google Moler predprey, it tries to talk me out of it, 167 00:13:20,250 --> 00:13:23,570 but then it shows it's in my second book, Experiments 168 00:13:23,570 --> 00:13:24,500 with MATLAB. 169 00:13:24,500 --> 00:13:26,000 There are two books. 170 00:13:26,000 --> 00:13:30,410 Numerical Computing with MATLAB and Experiments with MATLAB. 171 00:13:30,410 --> 00:13:35,190 And pred prey is in the second one, Experiments with MATLAB. 172 00:13:35,190 --> 00:13:37,640 You can go to website and you can 173 00:13:37,640 --> 00:13:42,890 download all of the programs from EXM 174 00:13:42,890 --> 00:13:50,830 or you can go down here and here's pred prey. 175 00:13:50,830 --> 00:13:59,001 So it's in Experiments with MATLAB on the MathWorks 176 00:13:59,001 --> 00:13:59,500 website. 177 00:14:05,390 --> 00:14:10,980 Just Google Moler predator prey and you'll find it.