cs222 »

**These are draft notes from subtitles. Please help to improve them. Thank you!**

Contents

- 1 Unit 4: Responsible Fishing
- 1.1 01 l Welcome Unit 4
- 1.2 cs222 unit4 01 q Simple Harvesting2
- 1.3 cs222 unit4 01 s Simple Harvesting2
- 1.4 cs222 unit4 02 q Logistic Growth2
- 1.5 cs222 unit4 02 s Logistic Growth2
- 1.6 cs222 unit4 03 p Fish Modeling2
- 1.7 cs222 unit4 03 s Fish Modeling2
- 1.8 cs222 unit4 04 l Equilibrium2
- 1.9 cs222 unit4 05 q Finding Equilibria2
- 1.10 cs222 unit4 05 s Finding Equilibria2
- 1.11 cs222 unit4 06 q Classifying Stability2
- 1.12 cs222 unit4 06 s Classifying Stability2
- 1.13 cs222 unit4 07 q Changing Harvest Rate2
- 1.14 cs222 unit4 07 s Changing Harvest Rate2
- 1.15 cs222 unit4 08 p Ramping Up2
- 1.16 cs222 unit4 08 s Ramping Up2
- 1.17 cs222 unit4 09 l Slope Fields2
- 1.18 cs222 unit4 10 p Infinity Time2
- 1.19 cs222 unit4 10 s Infinity Time2
- 1.20 cs222 unit4 11 p Another Unless2
- 1.21 cs222 unit4 11 s Another Unless2
- 1.22 cs222 unit4 12 l Optimization2
- 1.23 cs222 unit4 13 p Linear Ramp2
- 1.24 cs222 unit4 13 s Linear Ramp2
- 1.25 cs222 unit4 14 l Sources of Error2
- 1.26 cs222 unit4 15 p Most Critical Parameter2
- 1.27 cs222 unit4 15 s Most Critical Parameter2
- 1.28 cs222 unit4 16 l Model Problems2
- 1.29 cs222 unit4 17 q Two Types of Fish2
- 1.30 cs222 unit4 17 s Two Types of Fish2
- 1.31 cs222 unit4 18 l Conclusion2
- 1.32 cs222 hw4 01 p By Catch Extinction
- 1.33 cs222 hw4 01 s By Catch Extinction
- 1.34 CS222_hw4_02_p_Food
- 1.35 cs222 hw4 02 s Food Chain

We're looking at two different processes that counteract each other. The growth of the fish and the harvesting of the fish. Let's ignore harvesting for the moment. A simple assumption about the growth would be that the amount of fish increases year by year by the same percentage similar to compound interest. To make things simple, we are measuring the amount of fish as a number of tons not by counting. Here is a quiz for you. Let's not include harvesting. What should the right of change of the amount of fish be: 0.5 per ton times the amount of fish minus 4 * 10⁵ years per ton, or minus 0.5 per year times the amount of fish plus 4 * 10⁵ tons per year, or plus 0.5 per year times the amount of fish minus 4 * 10⁵ tons per year.

The first answer can't be true because the units of measurement are wrong. This would be tons divided by tons which would leave us with a measurement unit of 1 and years per ton won't work either--we had to have tons per year, the rate of change of a number of tons, so this doesn't work. This expression would not more the growth but the decay of fish minus the constant times the number of fish and this would be something like immigration not harvesting. Some constant is being edit to the rate of change. We have an immigration of 4 times 10⁵ tons per year, so these two can't be the answer, and we are left with the third one--0.5 per year times the number of fish. The positive growth rate as is appropriate for exponential growth and then there is a loss with a minus sign, which models harvesting.

Of course, exponential growth is highly unrealistic. Let's have a look at this diagram. Now the horizontal axis is not the time, it's the amount of fish, and the vertical axis is the growth rate. With exponential growth, the growth rate does not depend on the amount of fish. Even if the ocean were full of fish, we would still have the same growth rate, which cannot be true. There should be some sort of carrying capacity at which the growth rate becomes zero. The easiest way to include this effect is to use the growth rate. that depends on the amount of fish in the linear fashion. This is called logistic growth. If the amount of fish is very small, we start with the high growth rate and as the amount of fish increases, the growth rate drops and if the amount of fish reaches the carrying capacity of the environment, there won't be any growth anymore. How do we turn this effect plus constant harvesting into an equation. Here are three options for you to choose from. Should the rate of change of the amount of fish be of 0.5/1 year times the amount minus 10⁶ tons minus 10⁵ tons per year or should it be 0.5 per year times 1 minus the amount over 10⁶ tons times the amount minus 10⁵ tons per year or should it be 0.5 per year times the amount minus 10⁶ tons times the amount minus 10⁵ tons per year.

With the first one, the units of measurement would be okay for a rate of change of tons per year. If you look at this first expression, you see that it's not doing what it supposed to be doing. Actually we are still multiplying the current amount of fish by 0.5 per year but still the same growth rate all of the time. The second term 0.5 per year times minus 10⁶ tons is just similar to harvesting a huge the amount of harvesting--this cannot be the solution. The third option doesn't work because the unit of measurement are wrong. Tons times tons divided by year, tons squared per year. That's not what we want to have. The second one is the one that works. The initial rate of change that's 0.5 per year is multiplied by a factor that's 1 if the amount of fish is 0--1 minus 0 over something and is multiplied by 0 if that amount is 10⁶ tons 1 minus 10⁶ tons/10⁶ tons. This factor modifies the initial rate of change in the way that we need it and here of course, we don't have harvesting.

Now let's use the computer to learn what happens with different harvest rates. If we would not be harvesting at all, we would expect to see something like an exponential growth at the beginning, but after a while the amount of fish would approach the maximum carrying capacity. If you were harvesting too strongly, you would expect the number of fish to quickly go down to zero and then of course stay zero--that's clearly over fishing. Now, implement this using the forward Euler method and make sure that you stop harvesting when the amount of fish reaches zero. All constants that you need are already set up in the code. You only have to provide the implementation of the forward Euler method to fill the array called fish.

Let's first have a look at the output. You see that with rates of 20,000 and 50,000 tons per year, we do not seem to be over fishing. However, with rates of 100,000 and 200,000 tons per year, we clearly are over fishing, fish go extinct. The line of code that actually implements the forward Euler method is produced straightforward. The harder part is to stop harvesting and this solution is as follows-- we check whether or not the amount of fish for the next step is equal to zero or negative. If so, we set a flag that tells us whether I should be using zero for the amount of fishing in all upcoming steps. Note that you cannot simply check whether the amount of fish for the next step is exactly equal to zero--that's probably not going to happen.

And these simulations for harvesting fish, you see several situations in which the amount of fish stays constant. You could say the falls up and the falls down exactly balance each other. This has caused an equilibrium or steady state or fixed point. These three terms are almost interchangeable. I'm going to use the term equilibrium. In most cases, equilibria can be classified as either stable or unstable even though in some cases you may have a mixed form. This position at the bottom of the valley is a stable equilibrium. If you start close to it, you will end up at that position. This position on top of the hill is an unstable equilibrium. If you happen to be in this precised spot, you are going to stay there, but if you're off just by a little, you're going to roll down the hill.

Let's look at this equation for just the growth and constant harvesting and let's take care of the fact that we cannot harvest any longer once the fish become extinct. At which amount of fish do we find an equilibrium? All of these are given as multiples of 1 million tons. 0, 0.33, 0.50, 0.55, 1.0, 1.21, 1.45 and 2.0.

Obviously, there is an equilibrium at 0--f stays 0 once it becomes 0. This is not nothing but an equilibrium and there's two more-- one is at 0.55 million tons and the other is at 1.45 million tons. The first one is obvious, but now let's look into where these two other equilibria come from. If we reach an equilibrium, the amount has to become constant. The rate of change becomes zero, which means that the gain per time has to balance the loss per time. The loss per time is constant but the gain per time depends on the current amount. Let's build this type of diagram to show how the gain depends on the amount. If you look closely, you see that this has to be a parabola. If the amount of fish is zero, we are multiplying by zero here and the gain will be zero. That's no surprise. If there is no fish, there is no growth. We know this point. When the amount of fish is equal to the maximum carrying capacity, this factor becomes zero. The product is zero again so we know this point. As a parabola is symmetric, the maximum has to occur in the middle between the zeros. This maximum sits at an amount of 1 * 10⁶ tons. We can plug this value 1 million tons into this expression and get that the gain here is 2.5 * 10⁵ tons per year, which is a little more than the loss, which is 2 * 10⁵ tons per year. We have one equilibrium here and one equilibrium there. At both of these points the gain and the loss, the gain and the loss balance each other. If you want to, you can solve this quadratic equation to find that this equilibrium sits at 0.55 * 1 million tons and this equilibrium sits at 1.45 * 1 million tons.

Now we know that there are three equilibria--one trivial, one at 0, and 2 further ones down the road. Which of these are stable and which of these are unstable?

The one at zero is stable, the next one is unstable, and the last one is stable again. Why is that the case? Assume that we are close to zero so that we have a small amount of fish and this harvesting rate is going to quickly reduce that amount to zero again, which makes zero stable. By the way, we can't approach zero from the other side. There won't be any negative amount of fish. The second equilibrium is unstable. Imagine that we had a little less fish that we would harvest more than we gain in this specific period, which means we are going to further drive down the number. We are moving away from that equilibrium. Similarly if you have a little more fish then the growth rate exceeds the harvesting rate. We are going to have even more fish. In both directions, we are moving away from the equilibrium It's unstable. And now the third equilibrium, if we have a little less fish, the growth rate being curve exceeds the harvest rate, more growth than harvest. We are going to end up with even more fish. We're moving up on the horizontal axis towards the equilibrium. If we have a little more fish, the harvest rate exceeds the growth rate. We are fishing more than we gain, which curves down the amount. So from both sides, we are going to approach this equilibrium, which makes it stable. If we are in this region, we are definitely over fishing. The amount of fish is going down all of the time towards zero. The harvesting rate is way to large for the small amount of fish that we have in this region.

Now we look at different harvest rates--2.0, 2.5, 3.0 times 10⁵ tons per year. The middle one happens to be the same 2.5 * 10⁵ tons per year. The question is will the amount of fish eventually decay to 0. Do we have a strong issue with over fishing? This is the quiz for you--what about the amount of fish for the blue harvest rate-- will it never shrink to 0, does it depend on our circumstances whether or not it shrinks to 0, will it almost surely decrease to 0, or will it even surely decrease to 0. The same for the hour's rate of 2.5 * 10⁵ tons per year and the same for the harvest rate of 2.0 * 10⁵ tons per year. This is going to lead us to the concept of maximum sustainable yield.

We have already covered the lowest harvest rate of 2.0 * 10⁵ tons per year, and we saw that we have a problem with over fishing if we are below this equilibrium point. If we are above that equilibrium point, the amount of fish will not decrease to zero. For the lower value, it depends. Namely, it depends on the initial amount of fish. The upper harvest rate clearly amounts to over fishing independent of the initial amount of fish. The loss is always larger than the gain. In this case, we surely have the situation that the amount of fish decreases to zero. If we use this harvest rate that just touches the peak of our gain graph, things get a little tricky. If we start with an amount of fish that's below 1 million tons, we lose more than we gain per time, which means that the amount of fish is going to become zero eventually If we are in the region above 1million tons also the loss is larger than the gain so we are moving towards lower amounts, but in principle, we would be stopping at this equilibrium point where the loss equals the gain. This is what would be called semi-stable equilibrium. On the one side, it is unstable and on the other side it's stable, but imagine what happens if we have eventually reached this equilibrium from the other side, just a tiny disturbance would suffice to make you drop down to the left end of that graph. Almost surely the amount of fish will decrease to zero if we pick the rate of 2.5*10⁵ tons per year. And this rate is called the maximum sustainable yield. In theory, you could be harvesting at this rate forever and ever, but of course, that's pretty dangerous whenever there's this largest disturbance.

We saw that fishing with the harvest rate that's equal to the maximum sustainable yield is pretty dangerous. Let's try to come up with something that's safer-- for instance, ramping up the harvest rate over the course of time. We don't harvest during the first four years and after a time of 6 years, we're not using the maximum sustainable yield but we're using something that's a little smaller 0.8 times the maximum sustainable yield. The 0.8 would be a safety factor and between 4 years and 6 years, let's have a linear ramp. Now implement this ramp in the computer model. The forward Euler's method is already in place. You just have to provide the code that computes the harvest rate, which is used here.

We can solve this for instance by distinguishing between three different cases. First case, time is between 0 and 4 years, the output is 0. Second case, the time is 6 years or greater. The output will be 0.8, the maximum sustainable yield. Third case, we need a function to describe this line if the time is between 4 years and 6 years. The slope of this line is 0.8 times the maximum sustainable yield divided by this time span 6 years minus 4 years. The slope is multiplied by t, but if we're just using this expression, this would be the result that we get. A line would intercept zero. The trick is to subtract 4 years from that time, which would make it work. We still have this expression still describes the line with the right slope, but if you plug in 4 years, this coefficient times 0 and thus the result is 0 for 4 years. In the resulting diagram, you see how the growth of the population softly bends down between 4 and 6 years,but we are not over fishing. The implementation makes use of a factor to modify the maximum harvest rate which is 0.8 times the maximum sustainable yield. If the time is larger than 6 years, we are using the factor of 1. If the time is not larger than 6 ix years but it's larger than 4 years, we're using the linear ramp, and otherwise, we are using the factor of 0.

This type of diagram--called a slope field--helps to understand why everything is determined once we fix a certain initial amount of fish in this case. The general term would be initial value. These arrows point in the direction in which our solution curve has to run at every point in time and every amount fish we could draw such an arrow. This arrow specifies the rate of change. Here the rate of change is slightly positive. Here it's approximately zero. Here it's slightly negative and here it is really hugely negative. To see what this system is doing, we simply follow these arrows or take this one here as an initial value. We first have to go down and then move from between here. Once you have specified the initial condition to initial value in this case the initial amount of fish, you can trace out such a solution curve that describes the dependency of the amount of fish on time or whatever your differential equation is about. There is such a solution curve and there is no second solution curve. This would not work. This is what mathematicians mean when they speak of existence and uniqueness of the solution. Once you specified the initial condition, there is a solution curve and there is only one solution curve. Well, unless I have to say first thing that may happen is that the solution escapes infinity in finite time then we don't know how to continue afterward or the solution may fall into a hole of the definition set of our equations. For instance imagine that we are simulating a machine and we do not have any idea about what's happening to the machine in this region then we can't continue either. The following quiz shows that this escape to infinity in finite time can happen easily.

This looks like a pretty harmless differential equation. We want the derivative of x with respect to time to be exrep as 1/2 and the initial value-- you know we have to specify the initial value to get one single solution and the initial value should be 3. Implement the forward Euler's method for that differential equation and then experiment what setting the ends time of that simulation, try to find out what the time of explosion of that solution is and use that in here with two decimal places as you submit the solution.

The Euler forward method is pretty simple again and the time that you should find for the explosion should be 0.64. This is a diagram resulting from that setting and this really makes it obvious that the solution is escaping to infinity and a finite amount of time. Actually, this behavior is not that much of a surprise. X always has to grow because the right hand side is always positive. The square is always positive and if you're adding to it and dividing by a positive number, but as x grows, the square is going to become really huge, so the rate of change becomes really huge, which leads to x become even higher and so and so on. This is going to escape to infinity in a finite amount of time. Footnote for the mathematically proficient, what we're getting here is actually the shifted version of the tangent function.

To finish this discussion on existence and uniqueness of solutions of differential equations, here is a final item from the mathematical cabinet of curiosities another unless, if you will, the derivative of x with respect to time should be the cube root of x, and if we're going to try out three different initial conditions--0, 0.01, and a really, really tiny number 10^-300. Implement the forward Euler's method for this differential equation and then have a look at what happens with these three different initial values.

Even though the overall equation is simple, there's two things to be careful about. First, you cannot write 1/3 because that would be 0. Integer divided by integer would be an integer so you have to put in at least one decimal point here. The second thing to be careful about is Python's power operator. If you forget the parenthesis here, it would be forming the first power and would be dividing the result by 3, which is not what we want obviously, and here's the result. If you start with an initial value of 0, the results stays at 0 all of the time, but if we start with a positive number that's ever so slightly above 0, such as 10⁻³⁰⁰, we get a completely different evolution and as you can see, there's hardly any difference between starting with 10⁻³⁰⁰ and starting with 0.01. From a numerical view point, this differential equation looks as though 0 is a highly unstable fixed point. Once you get a tiny bit above 0, the solution is going to quickly divert from that and mathematically speaking, this is a differential equation whose solution curves are not unique. If you start with 0, you will branch off at any point in time in this fashion and get another solution.

Now let's return to fishing but look at it from a different perspective-- the perspective of optimization. We could run our simulation with different but constant harvest rates and check what the total amount of harvest in 10 years would be. If we are using to a high a harvest rate, we are over fishing and hence lose from the total amount that would be possible. If we harvest too little, we are under fishing which is also not optimal. Optimization aims at finding the right value for the harvest rate-- the value that would maximize a total amount of harvest. In technical terms, the function that we are optimizing is called the objective function, and the harvest rate is a parameter it depends on. The optimization should find the best value of that harvest rate, meaning the value for which the objective function is maximum. The same would have to happen for instance if we look at the numbers of units sold in some simulation or the money earned in some simulation. Most probably, we want to maximize that value, but there's also a large number of optimization problems in which we want to minimize the objective function. For instance, we particularly want to minimize the duration of a trip or the money spent.

Finding the best value for a single parameter is a pretty easy thing to be doing. Let's try something more difficult--optimizing the linear ramp, that is finding the best value for the point of time at we which we start ramping out the value. the start may be earlier, the end may be earlier, the end maybe later, the start maybe later, and so on. We want to find out which of this curves would produce the highest output. Not every combination of these two parameters make sense. We do not want this start time to be larger than the end time. That won't really work, so in this diagram the combinations below the diagram now are forbidden, and what our program will be doing, it will be varying both parameters between 0 and 10 years only in the allowed region and it will put dots in that region that show the total amount of harvest for all of these combinations of parameters. this blue dot would mean that this specific combination of parameters leads to a rather large total amount, and this smaller dot would mean that this combination of parameters leads to a smaller total amount and so on. The stepping of these two parameters that control the ramp process is already done. Your job is to compute the total amount of harvest and return it in this variable.

To compute the total harvest is simply take care of by how much it changes in the current step then add this to the total amount. This is the result and you can see that there is not a single combination of start and end time, but rather a line of combinations, all of which work similarly well. So if you transform this 2-dimensional view into a 3-dimensional view, it would look like the reach of a mountain, not a single peak of a mountain.

To apply stimulation in real life and to apply optimization in real life, we have to have an idea about the sources and the amounts of errors that we are introducing. The model as such may not be good enough. For instance, we may be missing some vital effects in that model. The values that we are using for the model's parameters may be wrong. We may be using the wrong growth rate or we may be using the growth rate that's too big or we may be using a maximum carrying capacity that's too large. The initial data may be off. Do we really know the amount of fish we're starting with? And we may have numerical issues, maybe due to round off errors or maybe, which is more probable, due to finite step sizes. We look at numerics in Units 2 and 3. I'll discuss the problems with this model at the end of these units. What I want to consider now is how to study errors concerning the parameters and initial data. The nice thing about these is that we can check them in the model. We can use a computer to get an idea about the effects at least of errors we are introducing here, and the computer can help us to find out which of these parameters and initial data is the most critical concerning its range of uncertainty. If we try to get more certain data, which of these would we try to improve? That's one of the major objectives of sensitivity analysis.

The most straightforward way to do sensitivity analysis is to do it one at a time. One factor at a time, that is. For instance, there is an uncertainty concerning the initial data. Actually, we may start here or we may start there and our knowledge is to where precisely we start may be limited. There is an uncertainty concerning the initial data, or we may not be able to precisely control the start of the ramp. If we start a little early, the total number of fish will be lower, as if we start a little late. As we change the single factors, sometimes the result goes up. As we change a single factor for instance, as in this case with the initial value, and sometimes as we look at the upper limit of one of the factors, the result may go down and it may go up as we look at the lower limit. Like in this case for the time at which the ramp starts. We can apply this to all factors that influence the result and check the degree by which they influence the result. I am just using the base value and then for every factor changing that factor to its upper limit and then to its lower limit. That's a little lazy. At times, this can be too lazy. Most often, the other possible values for that factor will lead to sensible curves in that range, but it may happen that sometimes you get very weird results for values in between. Nonetheless, in the upcoming program, we're looking at the outermost values. This code makes heavy use of dictionaries, a special feature of the Python language. A dictionary associates keys and values, so I can ask for the string initial_value and get 2 * 10^e5. I can ask for the string maximum_growth_rate and retrieve the value 0.5. This is a nice trick to have all factors in one place, because we need to change them one at a time. These are the base values and here we specify the uncertainties. So what we're saying here is that we know that the initial value of the amount of fish is between 2 * 10⁵ - 5 * 10⁴ tons and 2 * 10⁵ + 5 * 10⁴ tons and so on and so on. This dictionary is simply to make things colorful. The different factors will get different colors as well in the diagram. If you look at the code that we provide, you'll see more details about dictionaries in there. In other programming languages, you may find dictionaries too, or you may have to look for associative arrays or for maps or for hash tables. From the information provided in the code, you should be able to return the key of the most critical parameter in this variable.

The diagram that we get works as follows. Let's for instance look at the two red curves. They show the effect the initial value has on the result. It changes the initial value up and down by 5 * 10⁴ tons and you get this curve and that curve. You can immediately spot that this blue curve belongs to a factor that's more critical. The blue one is the carrying capacity, and obviously, the most critical is M for magenta, the maximum growth rate. If we want to improve the reliability of our result, we definitely need to work on the uncertainty of all maximum growth rate as opposed to what we just did. This code doesn't look at the complete temporal evolution, it just looks at one single value, the total amount of harvest for that particular set of parameters. It compares what happens when that factor goes up to what happens when that factor goes down and then determines the factor for which this difference is largest. When you show here as we increase the value of that factor, the result may go down, and as we decrease the value of that factor, the result may go up, so that we end up with a negative sign of the difference here. This is why we're using the absolute value and corresponding to what we guessed from the curves, the most critical parameter is the maximum growth rate.

Let's take a second look at the sources of errors. Mathematics can help us find and solve numerical issues. We can use the computer to a evaluate our model concerning uncertainty and sensitivity but what about the model as such? The model that we are using is really myopic. It completely ignores external influences. It doesn't cover spatial variation. What's going to happen if I over fish these fishes but I'll leave these fishes alone for instance? It doesn't include by-catch or if I try to harvest one species of fish, I will almost invariably also find other fishes in the net including some of that effect is one of homeworks, but the most drastic omission is that we don't think about food webs. Consider for instance this situation, the purple species of fish feeds on plankton but at the same time, these purple species is prey to the blue species. If you harvest the blue species of fish, the purple species of fish will grow in volume and may deplete the plankton and of course, there is a lot more of issues with this model.

Let's improve our myopic model a little. What could these two differential equations stand for that fish of type 1 die, there are young fish and mature fish, and that fish of type 2 are harvested. For each of these options, pick yes, maybe or no.

Welcome to the first problem of Unit 4. In this situation, we're going to pretend that we are fishermen. And we're interested in catching this red type of fish right here. However, whenever we try to catch the red fish we also end up catching a bunch of these smaller green fish, as you can see from the contents of the net right here. What we're going to do in this problem is model logistic growth for both species of fish-- the red and the green with the constants and initial values that we've already included in the supplied code. Just a reference--I included a general equation from logistic growth right here. F is the population that we're looking at, r is the growth rate of the population, and c is the carrying capacity of the surrounding environment. Now, remember when we model the change in population for both types of fish, we need to include both population growth, which is dictated by this equation, but also the harvesting rate. We've created this parameter called p to represent the fraction of the total fish that we catch that are part of the green fish population, so we might also call this the fraction of the by-catch. Just to start, we've inserted a value of 0.5 for p. What you're going to do is calculate the threshold value of p for fish 2 to not go extinct. In other words, what is the highest possible value of p that will allow for equilibrium situation in the population of fish 2? After you've calculated the correct value of p, come down to the for loop and use, as we've done so many times before, the forward Euler method. You'll notice that we've already included a line for you to prevent the population of either type of fish from becoming negative--good luck.

Now let's look at the solution to the first problem. We know that we already used the forward Euler method for both the populations of fish one and fish two. We need to know how each population changes with time. I've written out expressions for ḟ₁ and ḟ₂ showing the time derivatives of each population. Now on each equation, we can see that this first term is in the form of the logistic growth function that I showed you in the question introduction video. We have a growth rate and then 1 minus the population over the carrying capacity times the population again and the same thing up here. These first terms show the positive rate of change of the population. However, each equation also has a second part tacked onto it that shows a negative rate of change. This expresses how each population is harvested by the fisherman. As I said earlier, p is the fraction of the fish that we catch that are fish type two or the green fish. We know then that the rate of harvesting of the green type of fish is going to be p times our total harvest rate. We can maintain the same total harvest rate if we make the coefficient in front of the harvest rate a new equation for the red fish 1-p. Remember we're trying to prevent the green fish from going extinct. This means that we need to create a situation in which the rate of change of population two is greater than or equal to zero. If we rearrange the equation for ḟ2, which should adapt with this term, which is showing how the population can gain fish, should be greater than or equal to this term, which is showing how the population can lose fish. The value of p that we're trying to calculate, however, is the exact threshold of extinction. This is going to happen where growth exactly equals harvesting. In other words if there are anymore harvesting or any less growth, population two will become extinct eventually. This graph over here like you've seen before has a blue pair of fish that shows how the growth rate changes in population, and it also has a red curve that represents harvesting rate. This first red line is just an example of a harvest rate, but the one that we're interested in is this dotted red line up here, which just touches the maximum of the growth curve. You might remember that this point on the growth curve is also known as the maximum sustainable yield. You can see that if we move the red line any higher, then the harvesting rate would always succeed the growth rate and we would definitely end up with an extinct fish two population. Similarly, if we move anywhere else on the blue line, the growth rate would again be below this harvest rate and we will begin to deplete the population of fish two. Right in the center, however, we have a point of semi-stable equilibrium just like we're going to talk about in the Unit. Any tiny bit that the population moves to the left side of the curve will result in it eventually falling all the way down to zero, which clearly means that fish two has become extinct. Instead, we have slightly larger population than 5x10⁵ tons, we will again start to see a decrease in population, and if by chance the population has not stopped at exactly this maximum point--it will fall to the other side of the curve and also become extinct. Remember all this is assuming that the harvest rate is exactly at this dotted line. Now we need to find what value of p would place the dotted line right at this position. Again at this point we have growth rate exactly equal to harvesting rate. I've taken the two parts of the equation for ḟ2 and written them as equal to each other right here. We plug in the proper value for f2, which we know is 5x10⁵ ton and we can solve using simple algebra for the threshold value of p, which equals to 0.375. Finally looking at our code, we can see that first we've plugged in the value for p that we just calculated, 0.375, and down here in the for loop, we have a very simple translation of the equations that we created for the growth rates of both fish one and fish two in this box for the rate of change of either population. So for fish one we have this and for fish two we have this. As always, these rates of change are multiplied by h, the step size, and added to the initial values for either fish population. Now it's our new program and see what we get. Here is our plots, the blue line plots the population fish that we are interested in harvesting measured in tons and the green line plots the population of the second type of fish, which we're trying to prevent from going extinct. You could see that we have maintained a high level of population of the fish that we're interested in and we've also prevented the green line from going down to zero. You can see that as time increases, the slope of this line is approaching zero, which means that if we plotted time going even farther out, no matter how far out in time we decided to to plot these curves, this green line would never approach zero. Just for fun, let's see what happens if we set p=0.5. Well this graph looks very different from the last one that we had. By just increasing p by the third of its former value, we've made fish population two go extinct in just over 20 years. Since we're harvesting more of fish two and less of fish one, the fish one population has actually risen to a value that is higher than it was in the last part, but we've violated our first main objective, which was to prevent fish two from going extinct. I hope this problem has helped you submit your understanding of logistic growth, harvesting, and equilibria.

Looking at the two model equations we gave you, you can figure out which one belongs to prey population and which one to the predator population by considering the signs of the various constants in the information that we gave you. We told you that the prey grows at constant rate of 0.5 per year. So since population growth is, of course, positive it depends on the population already present, we can see that this term Ax(t) should be part of the equation dealing with the prey. Similarly, we told you that the predators die out. Similarly, we told you that the predators have a fixed lifespan of five years. This means that the rate of change of the predator population while the negative component that is also directly proportional to the total population at a given time. This matches this term right here, -Cy(t). Now that we know that x represents predator and y represents prey, Let's talk about these terms that includes an x*y factor. This factor of x*y represents predator times prey and this proportional to basically how often individuals of the two species meet. We know that whenever a predator meets its prey, the situation, as it has been, pretty advantageous for the predator but not so great for the prey. This was shown in our equations by the fact that this interaction turn in the equation for the rate of change of the predator population has a positive sign in front of it. Whenever the predator meets its prey, this population can grow, and opposite is true for the prey population since this interaction will lead to members of the prey dying. Determining these constants A, B, C, and D should actually remind you of it of dealing with our SIR model problems in Unit 3. Our coefficients dealing with lifespan and growth rate corresponds to the different time constants that we use from moving people and mosquitoes to the infected, recovered, and susceptible populations. Since we've already matched with the growth rate of the prey population, we can write down A=0.5/year In the same way, since C corresponds to loss in the predator population, we can set its value equal to 1 over the predator lifespan or 1/5 years. In order to calculate B and D, the other two constants, we need to make use of the last piece of information that I gave you. That our ocean populations reach an equilibrium situation then there are 5.0*10⁶ tons of prey and 1.0*10⁶ tons of predator. We do some simple algebraic manipulations of the equations that I showed you above and also plug in 0 with the rates of change divided by population showing equilibrium situation. We can solve for B and D by simply inputting the values for A and C that we just decided upon. This results in B=5.0*10⁻⁷, 1/year*ton. D=4.0*10⁻⁸, 1/year*ton. We plugged all four of these values into base values. We can move on to the next part of the problem. In the function food chain, we implement the forward Euler method using the rates of change that we just figured out using our differential equations. Of course, you've plugged in prey in the place of x and predator in the place of y. And finally, we move on to our sensitivity analysis. To account for the lower limit values, we make a copy of the base values dictionary and multiply each value corresponding to each key by 0.9 and we do the same thing for the upper limit by multiplying each value by 1.1. These lines down here select the most critical parameter based on which parameter has the greatest difference between its high results and its low results. However, we can also figure out the most critical parameter by just looking at the plots that we get when we run the program. Our first plot shows the amount of prey in tons versus time and the second shows the amount of predator in tons versus time. The colors dictionary which you may have seen throughout the top of the code paired each of the code parameter that we're looking at with a different color. Remember, we're interested in seeing which of these parameters has the greatest impact on the maximum value of the prey? So in order to pick out the most critical parameter, we need to compare the maximum amount of prey graphed by the upper and lower curves of either series in this top graph. You can see for example that if we look at the yellow curves, their maximum values here and here are not really that far apart. If we look at every pair of curves that are the same color, we can see that this one that is cyan or the light blue has the greatest distance between its peaks. Looking back at our colors dictionary, we see that cyan or C is paired with D. Let's look at our final graph again for just a little bit more analysis. We can see that both graphs exhibit periodic behavior. It looks almost like we have a Cos and Sin function. Initially, our prey population increases which leads to an increase in the predator population shortly afterward. This however leads to a decrease in the prey population and so on and so forth back and forth. Upon closer inspection, you can see that these are actually not perfect Sin and Cos functions. In fact, the Euler method leads to a blow-up of periodic functions overtime. We see a very clear example of this explosion of amplitude if we increase the end time and step size. Let's see what our plot looks like now. This is a very interesting-looking plot. We can see that the in the predator population, we start with a peak-to-peak height of approximately 60,000 tons and after 200 years have passed, we have a peak-to-peak height of approximately 150,000 tons. That's a pretty dramatic increase. If we use this implicit Euler method instead of the forward Euler method here, we could prevent this expansion and amplitude just as we talked about in our earlier problems involving mechanics. You can see that many methods prove useful in different situations whether they're dealing with diseases or fish populations or pendulum. Great job on Unit 4 and get excited for Unit 5. We're going to look at something totally different. How anti-lock brake systems work.