Contents
Mathematical Model:...................................................................................................................................2
1. Classify (1) and interpret (in words) each term of f(y).........................................................................2
2. When interpreting H(y), consider these questions: what happens to H(y) as y gets very large or
very small? What does this represent physically? Does this make sense? (Hint: It may help to
remember that y0 and, therefore, H(y) has units of hundreds of tribbles per day.)...............................3
3. units of parameters a and b.............................................................................................................3
4. Solve analytically.................................................................................................................................4
5. Consider an ordinary differential equation (ODE) of the form y0 = f(y)...............................................7
6. For b = 0:005; 0:05; and 0:10, plot f as a function of y. Use this plot to determine (graphically) the
equilibrium solutions of (1) for each value of b. Make sure you use a suitable scale when plotting, so
that you do not miss any equilibrium values. You can also use a root finding routine to find the
equilibrium solutions (e.g. fzero in Matlab). Plotting f(y) in Desmos may also be beneficial (enter b as
a parameter/slider).................................................................................................................................8
7. Plot direction fields of (1) and overlay several solutions for the same values of b as in question 6.
Consider the equilibrium solutions when choosing your initial conditions for your solutions. You can
do this in MATLAB (quiver command will aid in creating the direction fields; you can code up Euler’s
method or use ode45 to find solutions) or you can use the Geogebra script located here which will
plot solutions and direction fields together. The files dirfield.m and basic ode45.m in the Canvas
Project Information module may also be of some help.........................................................................11
8. Interpret the function plots, direction fields and solutions.................................................................13
9. Suppose brown tribbles have b = 0.005, white have b = 0.05, and grey have b = 0.10. At the
beginning of tribbling season, we can ensure the initial population to be at a certain value. Based on
your results above, with which tribble would it be best to stock the station? How many tribbles should
we make sure are in the station at the beginning of tribbling season (using your previous species
choice)? Note: this question is intentionally vague and somewhat open-ended. Your job is to answer it
using the tools you have developed. You will be graded not only on how well you apply these tools, but
also on your clarity in stating priorities and assumptions and realizations (mathematically,
computationally, etc.) of those priorities and assumptions....................................................................13
10. Seasonal effect of harvesting...........................................................................................................16
11. Discuss briefly concerning the following issues. Are there any weaknesses in the model used? Are
there additional effects that the tribble population model should take into account? Can your model be
improved and if so, how? (These are just suggestions for some interesting questions. You don’t have to
address all of these questions in the report.............................................................................................19
Appendix...................................................................................................................................................20
Mathematical Model:
1. Classify (1) and interpret (in words) each term of f(y).
This term represents the contribution of the birth rate to the population growth. The parameter 'a'
determines the strength of the birth rate effect. The term ay indicates that the population size y is
multiplied by the birth rate constant 'a', resulting in an increase in the population size over time.
-by^2:
This term represents the contribution of the death rate to the population decline. The parameter 'b'
determines the strength of the death rate effect. The term -by^2 indicates that the population size y is
squared and multiplied by the negative death rate constant '-b', leading to a decrease in the population size
due to mortality.
-H(y):
This term represents the harvesting effect on the population dynamics. The function H(y) = py^3/(y^3 +
q) determines the harvesting rate, which depends on the population size y. The parameter 'p' determines
the strength of the harvesting effect, while 'q' affects the saturation level of the harvesting effect.
When interpreting H(y), consider these questions: what happens to H(y) as y gets very large or
very small? What does this represent physically? Does this make sense? (Hint: It may help to
remember that y0 and, therefore, H(y) has units of hundreds of tribbles per day.)
The harvesting term -H(y) implies that the harvesting rate is subtracted from the population
growth. As the population size increases, the harvesting effect becomes more pronounced due to
the cubic term py^3 in the numerator. The denominator y^3 + q ensures that the harvesting rate
reaches a saturation level as the population size grows, preventing the harvesting effect from
becoming unbounded.
As y gets very large, the term y^3 dominates the denominator y^3 + q, making the harvesting
term approach py^3/y^3. In this case, H(y) approaches p, meaning that the harvesting rate
becomes constant and independent of the population size. This suggests that when the population
of tribbles is large, the harvesting effect remains relatively consistent.
As y gets very small, the term q dominates the denominator y^3 + q. This causes the harvesting
term to approach py^3/q. Consequently, H(y) approaches zero as y approaches zero, indicating
that the harvesting effect becomes negligible for small population sizes.
2. When interpreting H(y), consider these questions: what happens
to H(y) as y gets very large or very small? What does this
represent physically? Does this make sense? (Hint: It may help to
remember that y0 and, therefore, H(y) has units of hundreds of
tribbles per day.)
Physically, the behavior of H(y) reflects the dynamics of the harvesting effect on the tribble
population. As the population size grows larger, the harvesting effect becomes more prominent
and reaches a saturation level determined by the parameter q. This saturation suggests that there
is a limit to the impact of harvesting, even as the population continues to increase.
Conversely, for small population sizes, the harvesting effect becomes less significant. This can
be interpreted as a reduced harvesting impact due to the limited availability of tribbles to harvest
or the decreased interest in hunting when the population is sparse.
Yes, the behavior of H(y) in relation to the population size y makes sense. As the population size
increases, the harvesting effect becomes more pronounced, which aligns with the idea that a
larger population provides more resources for harvesting. Similarly, as the population size
decreases, the harvesting effect diminishes, reflecting the reduced impact of harvesting on a
smaller population.
3. units of parameters a and b.
The left-hand side dy/dt represents the change in population per change in time, which has units of tribbles per day.
Similarly, we have
4. Solve analytically
a) Use separation of variables to solve (1) analytically when H(y) = 0, that is, when p = 0. Assume the initial
population y(0) = y0. What does a=b represent? (Hint: look at what happens to your solution as t grows without
bound)
As t grows without bound, it means birth rate is higher than death rate and in case of a= b , it means
population achieves steady state solution.
(b) If H(y) is non zero in (1), can we still solve it analytically? If so, how? Set up the integrals needed and explain
how they could be evaluated. Do not evaluate them
By partial fraction expansion, the polynomial term would be expanded and then integral will be
evaluated.
5. Consider an ordinary differential equation (ODE) of the form y0 = f(y)
(a) Define in words what an equilibrium solution to the equation is
An equilibrium solution to the equation y' = f(y) represents a state in which the dependent variable y
remains constant over time. In other words, at an equilibrium solution, the rate of change of y, represented
by y', is zero.
Mathematically, an equilibrium solution occurs when f(y) = 0. This means that the function f(y) does not
cause any change in the value of y, resulting in a steady state.
Interpreting this conceptually, an equilibrium solution can represent a stable point or state of balance in a
system. It is a solution where the forces or factors influencing y reach an equilibrium, resulting in a
constant value of y over time. At this point, the rate of change of y is zero, indicating that the system is
neither increasing nor decreasing.
(b) Give an equivalent mathematical definition and explain how you use this to find the equilibrium
solutions
Mathematically, an equilibrium solution occurs when y' = 0, which implies that the derivative of y with
respect to the independent variable is zero. This condition indicates that there is no rate of change in y,
and thus y remains constant.
To find the equilibrium solutions, we set y' = 0 and solve the equation f(y) = 0. This means that at an
equilibrium solution, the function f(y) evaluates to zero, indicating that the factors influencing the rate of
change of y have balanced out, resulting in a constant value of y.
To solve f(y) = 0, we can apply algebraic or numerical methods depending on the specific form of the
function f(y). These methods may include factoring, root-finding techniques, or solving transcendental
equations.
There can be multiple equilibrium solutions depending on the behavior of the function f(y).
6. For b = 0:005; 0:05; and 0:10, plot f as a function of y. Use this plot
to determine (graphically) the equilibrium solutions of (1) for each
value of b. Make sure you use a suitable scale when plotting, so that
you do not miss any equilibrium values. You can also use a root finding
routine to find the equilibrium solutions (e.g. fzero in Matlab). Plotting
f(y) in Desmos may also be beneficial (enter b as a parameter/slider)
Equatilibrium point is: 0, 147.973
Equatilibrium point is: 12.625
Equatilibrium point is: 1.077
Equatilibrium point is: 0
Equatilibrium point is: 0.000
Equatilibrium point is: 0.971
Equatilibrium point is: -1.077
7. Plot direction fields of (1) and overlay several solutions for the same
values of b as in question 6. Consider the equilibrium solutions when
choosing your initial conditions for your solutions. You can do this in
MATLAB (quiver command will aid in creating the direction fields; you
can code up Euler’s method or use ode45 to find solutions) or you can
use the Geogebra script located here which will plot solutions and
direction fields together. The files dirfield.m and basic ode45.m in the
Canvas Project Information module may also be of some help
Negative initial condition is just to show hypothetical equilibrium solution when initial condition attracts toward
zero trivial equilibrium solution.
8. Interpret the function plots, direction fields and solutions
Any initial condition above 0.01 leads to equilibrium value. The stable equilibrium value for each
of the different value of b is different.
With smallest b values related to least death rate means a larger equilibrium value. Like for
0.005, the equilibrium value is around 150. For a b value of 0.05 a larger death rate as compared
to 0.005, gives an equilibrium value of around 12 hundreds of tribbles. Similarly, a higher death
rate leads to and equilibrium value of unit hundreds of tribbles.
9. Suppose brown tribbles have b = 0.005, white have b = 0.05, and grey have b = 0.10.
At the beginning of tribbling season, we can ensure the initial population to be at a
certain value. Based on your results above, with which tribble would it be best to stock
the station? How many tribbles should we make sure are in the station at the beginning of
tribbling season (using your previous species choice)? Note: this question is intentionally
vague and somewhat open-ended. Your job is to answer it using the tools you have
developed. You will be graded not only on how well you apply these tools, but also on
your clarity in stating priorities and assumptions and realizations (mathematically,
computationally, etc.) of those priorities and assumptions
For this purpose, basin of attraction for each of the tribbles type that are brown, white and grey.
Brown tribbles
Equilibrium value is a/b = 0.75/0.005 = 150 hundred of tribbles
From above figure basin of attraction for brown tribbles is 1 to 150 and then infinity.
White tribbles:
Equilibrium value is a/b = 0.75/0.05 = 15.0 hundred of tribbles
From above figure basin of attraction for brown tribbles is 1 to 15 and then infinity.
Grey tribbles:
Equilibrium value is a/b = 0.75/0.05 = 1.5 hundreds of tribbles
From above figure basin of attraction for brown tribbles is 0.1 to 1 and then infinity.
Based on the calculations above, we can determine the best tribble species to stock the station by
considering the basin of attraction for each species. The basin of attraction represents the range of initial
population values that will eventually lead to a stable equilibrium population.
For brown tribbles, the basin of attraction ranges from 1 to 150 hundred tribbles. This means that if the
initial population is within this range, it will eventually stabilize around 150 hundred tribbles.
For white tribbles, the basin of attraction ranges from 1 to 15 hundred tribbles. Thus, if the initial
population is within this range, it will stabilize around 15 hundred tribbles.
For grey tribbles, the basin of attraction ranges from 1 to 1.5 hundred tribbles. So, if the initial population
is within this range, it will eventually stabilize around 1.5 hundred tribbles.
Considering these basin of attraction values, we can conclude the following:
If the goal is to maintain a large tribble population, the best choice would be brown tribbles, as they have
the highest equilibrium value and the widest range of initial population values that lead to stable
populations.
If the goal is to maintain a moderate-sized population, white tribbles could be a suitable choice, as they
have an intermediate equilibrium value and a narrower basin of attraction compared to brown tribbles.
If the goal is to maintain a small population, grey tribbles would be the most appropriate choice, as they
have the lowest equilibrium value and the narrowest basin of attraction.
10. Seasonal effect of harvesting
(b) What is the long-term behavior? What has changed in the problem? Can you analytically solve it? Are there any
equilibrium solutions left?
To analyze the long-term behavior of the given model: dy/dt = ay - by^2 - H(t, y), we need to consider the equation
without the time-dependent term H(t, y) first. The equation becomes a standard logistic differential equation:
dy/dt = ay - by^2
This logistic equation has two equilibrium solutions: y = 0 and y = a/b. The equilibrium solution y = 0 is stable if a >
0 and unstable if a < 0. The equilibrium solution y = a/b is stable if a > 0 and unstable if a < 0.
Now, let's consider the effect of the harvesting term H(t, y) = py^3 / (q + y^3) * abs(sin(pit/360)). This harvesting
term introduces a seasonal effect to the population growth. The term depends on both the population size (y) and
time (t). Analytically solving the modified model with the harvesting term is challenging due to the nonlinearity and
time-dependence of the term. However, we can still analyze the behavior qualitatively. The harvesting term H(t, y)
can have two effects: a positive effect when sin(pit/360) is positive and a negative effect when sin(pit/360) is
negative. The magnitude of the effect depends on the values of p, q, and the population size y. The positive effect of
the harvesting term can lead to an increase in the population growth rate, counteracting the negative effects of the
logistic growth term. Conversely, the negative effect of the harvesting term can reduce the population growth rate.
Considering the given data, we can observe that the population size (y) decreases over half of the cycle and then
start increasing after mid of the year. This suggests that the negative effects of the harvesting term dominate first
half and positive effects dominates in the next half of the year.
11. Discuss briefly concerning the following issues. Are there any weaknesses in the
model used? Are there additional effects that the tribble population model should take
into account? Can your model be improved and if so, how? (These are just suggestions
for some interesting questions. You don’t have to address all of these questions in the
report.
The tribble population model discussed earlier has some strengths but also several weaknesses. Let's discuss them
along with additional effects that could be considered, and potential improvements for the model:
Weaknesses in the Model:
1. Lack of Environmental Factors
The model assumes a constant growth rate without considering external factors such as resource availability,
predation, or competition. In reality, these factors can significantly impact population dynamics.
2. Limited Genetic Considerations
The model does not account for genetic factors, such as mutations or genetic variation, which can influence the
growth and adaptability of a population.
3. Homogeneous Population
The model assumes a uniform population, disregarding potential variations within the tribble population. Different
individuals may have different reproductive rates or survival capabilities, which could affect the overall population
growth. There is only one specie and one factor is death and birth there is no other factors included like hunting or
rivalry in the population.
Additional Effects to Consider:
1. Carrying Capacity
The model could incorporate the concept of carrying capacity, which represents the maximum population size an
environment can support. Beyond this threshold, the population growth rate would decline.
2. Age Structure
The age structure of the population can have a significant impact on growth rates. Incorporating different
reproductive rates for various age groups could provide a more accurate representation.
3. Interactions with Other Species
The model could consider the interactions between tribbles and other species, such as predators or competitors.
These interactions can influence population growth and survival.
Improvements to the Model:
1. Environmental Feedback
Introducing feedback mechanisms that consider resource availability, predation, and competition would enhance the
model's realism. For example, the growth rate could be adjusted based on the availability of food or the presence of
predators.
2. Interactions with Other Species
Like predator prey model, some specie will be required as prey and some acts as predator in the food chain.
Considering it will definitely improve the model.
Nonetheless, addressing these weaknesses and incorporating additional effects would result in a more realistic and
robust tribble population model.
Appendix
Project 1 : Tribbles population
clear all
close all
clc
%For b = 0:005; 0:05; and 0:10, plot f as a function of y. Use this plot to determine (graphically) the equilibrium
solutions of (1) for
%each value of b. Make sure you use a suitable scale when plotting, so that you do not miss any equilibrium values.
You can also use a
%root finding routine to find the equilibrium solutions (e.g. fzero in Matlab). Plotting f(y) in Desmos may also be
beneficial (enter b
%as a parameter/slider)
a =0.75;
bs = [0.005,0.05,0.1];
b = bs(1);
p = 1.5;
q = 1.25;
f = @(y) a*y -b*y.^2 -p*y.^3/(q+y.^3)
figure(1)
hold on
fplot(f)
hold on
%plot(y,fys)
plot(-1:1:200,zeros(length(-1:1:200),1))
xlabel('x')
ylabel('f(y)')
title('f(y)=ay-by^2-H(y) with a = 0.75 , b = 0.005, p = 1.5, q = 1.25')
grid on
f=
function_handle with value:
@(y)a*y-b*y.^2-p*y.^3/(q+y.^3)
Warning: Function behaves unexpectedly on array inputs. To improve performance,
properly vectorize your function to return an output with the same size and
shape as the input arguments.
finding roots
yeq=fzero(f,100);
fprintf("Equatilibrium point is: %.3f \n",yeq)
Equatilibrium point is: 147.973
b = bs(2);
p = 1.5;
q = 1.25;
f = @(y) a*y -b*y.^2 -p*y.^3/(q+y.^3)
figure(2)
hold on
fplot(f)
hold on
%plot(y,fys)
plot(-1:1:200,zeros(length(-1:1:200),1))
xlim([-5 25])
xlabel('x')
ylabel('f(y)')
title('f(y)=ay-by^2-H(y) with a = 0.75 , b = 0.005, p = 1.5, q = 1.25')
grid on
f=
function_handle with value:
@(y)a*y-b*y.^2-p*y.^3/(q+y.^3)
Warning: Function behaves unexpectedly on array inputs. To improve performance,
properly vectorize your function to return an output with the same size and
shape as the input arguments.
finding roots
yeq=fzero(f,100);
fprintf("Equatilibrium point is: %.3f \n",yeq)
yeq=fzero(f,1);
fprintf("Equatilibrium point is: %.3f \n",yeq)
Equatilibrium point is: 12.625
Equatilibrium point is: 1.077
b = bs(3);
p = 1.5;
q = 1.25;
f = @(y) a*y -b*y.^2 -p*y.^3/(q+y.^3)
figure(3)
hold on
fplot(f)
hold on
%plot(y,fys)
plot(-1:1:200,zeros(length(-1:1:200),1))
xlim([-5 25])
xlabel('x')
ylabel('f(y)')
title('f(y)=ay-by^2-H(y) with a = 0.75 , b = 0.005, p = 1.5, q = 1.25')
grid on
f=
function_handle with value:
@(y)a*y-b*y.^2-p*y.^3/(q+y.^3)
Warning: Function behaves unexpectedly on array inputs. To improve performance,
properly vectorize your function to return an output with the same size and
shape as the input arguments.
finding roots
yeq=fzero(f,0);
fprintf("Equatilibrium point is: %.3f \n",yeq)
yeq=fzero(f,1);
fprintf("Equatilibrium point is: %.3f \n",yeq)
yeq=fzero(f,-1);
fprintf("Equatilibrium point is: %.3f \n",yeq)
Equatilibrium point is: 0.000
Equatilibrium point is: 0.971
Equatilibrium point is: -1.077
Part 7 directional field
a= 0.75;
b = 0.005;
p = 1.5;
q = 1.25;
% Define the range of y and t values
y_min = -2;
y_max = 200;
t_min = 0;
t_max = 100;
% Define the number of grid points
num_points = 20;
% Generate grid points
yval = linspace(y_min, y_max, num_points);
tval = linspace(t_min, t_max, num_points);
yprime=@(t,y) 0.75.*y -0.005.*y.^2 -1.5.*y.^3./(1.25+y.^3);
figure(4)
%dirfield(yprime,tval,yval)
hold on
dirfield(yprime,tval,yval)
y0=1;
[t y]=ode45(yprime,[0 100],y0);
plot(t,y,'-k','DisplayName',['y_0 = ',num2str(y0)],'linewidth',2)
y0=0.1;
[t y]=ode45(yprime,[0 100],y0);
plot(t,y,'-b','DisplayName',['y_0 = ',num2str(y0)],'linewidth',2)
y0=10;
[t y]=ode45(yprime,[0 100],y0);
plot(t,y,'-m','DisplayName',['y_0 = ',num2str(y0)],'linewidth',2)
y0=200;
[t y]=ode45(yprime,[0 100],y0);
plot(t,y,'-r','DisplayName',['y_0 = ',num2str(y0)],'linewidth',2)
y0=-1;
[t y]=ode45(yprime,[0 100],y0);
plot(t,y,'-m','DisplayName',['y_0 = ',num2str(y0)],'linewidth',2)
title('Direction fields plot for a = 0.75, b = 0.005, p = 1.5, q=1.25')
legend('location','northeast')
b = 0.05
yprime=@(t,y) 0.75.*y -0.05.*y.^2 -1.5.*y.^3./(1.25+y.^3);
% Define the range of y and t values
y_min = -2;
y_max = 20;
t_min = 0;
t_max = 100;
% Define the number of grid points
num_points = 20;
% Generate grid points
yval = linspace(y_min, y_max, num_points);
tval = linspace(t_min, t_max, num_points);
figure(5)
%dirfield(yprime,tval,yval)
hold on
dirfield(yprime,tval,yval)
y0=1;
[t y]=ode45(yprime,[0 100],y0);
plot(t,y,'-k','DisplayName',['y_0 = ',num2str(y0)],'linewidth',2)
y0=0.1;
[t y]=ode45(yprime,[0 100],y0);
plot(t,y,'-b','DisplayName',['y_0 = ',num2str(y0)],'linewidth',2)
y0=10;
[t y]=ode45(yprime,[0 100],y0);
plot(t,y,'-m','DisplayName',['y_0 = ',num2str(y0)],'linewidth',2)
y0=200;
[t y]=ode45(yprime,[0 100],y0);
plot(t,y,'-r','DisplayName',['y_0 = ',num2str(y0)],'linewidth',2)
y0=-1;
[t y]=ode45(yprime,[0 100],y0);
plot(t,y,'-m','DisplayName',['y_0 = ',num2str(y0)],'linewidth',2)
title('Direction fields plot for a = 0.75, b = 0.05, p = 1.5, q=1.25')
legend('location','northeast')
b = 0.1
yprime=@(t,y) 0.75.*y -0.1.*y.^2 -1.5.*y.^3./(1.25+y.^3);
% Define the range of y and t values
y_min = -2;
y_max = 5;
t_min = 0;
t_max = 100;
% Define the number of grid points
num_points = 20;
% Generate grid points
yval = linspace(y_min, y_max, num_points);
tval = linspace(t_min, t_max, num_points);
figure(6)
%dirfield(yprime,tval,yval)
hold on
dirfield(yprime,tval,yval)
y0=1;
[t y]=ode45(yprime,[0 100],y0);
plot(t,y,'-k','DisplayName',['y_0 = ',num2str(y0)],'linewidth',2)
y0=0.1;
[t y]=ode45(yprime,[0 100],y0);
plot(t,y,'-b','DisplayName',['y_0 = ',num2str(y0)],'linewidth',2)
y0=10;
[t y]=ode45(yprime,[0 100],y0);
plot(t,y,'-m','DisplayName',['y_0 = ',num2str(y0)],'linewidth',2)
y0=200;
[t y]=ode45(yprime,[0 100],y0);
plot(t,y,'-r','DisplayName',['y_0 = ',num2str(y0)],'linewidth',2)
y0=-1;
[t y]=ode45(yprime,[0 100],y0);
plot(t,y,'-m','DisplayName',['y_0 = ',num2str(y0)],'linewidth',2)
title('Direction fields plot for a = 0.75, b = 0.1, p = 1.5, q=1.25')
legend('location','northeast')
function [yp] = rhs(t,y) % do not change this line
% *** Your entry 6 *** (start entering your system here)
a= 0.75;
b = 0.005;
p = 1.5;
q = 1.25;
yp=a*y -b*y.^2 -p*y.^3./(q+y.^3);
end
Project 1 : Part 10....................................................................................................................................1
b = 0.05....................................................................................................................................................2
Project 1 : Part 10
clear all
close all
clc
tspan=[0 365*10];
y0s= [0.01,0.1,1,10,100,150,250];
figure(1)
yprime=@(t,y) 0.75.*y -0.005.*y.^2-1.5*y.^3./(1.25+y.^3).*abs(sin(pi.*t/365));
% Define the range of y and t values
y_min = 0;
y_max = 200;
t_min = 0;
t_max = 365*10;
% Define the number of grid points
num_points = 100;
% Generate grid points
yval = linspace(y_min, y_max, num_points);
tval = linspace(t_min, t_max, num_points);
hold on
dirfield(yprime,tval,yval)
for i=1:1:length(y0s)
y0=y0s(i);
[t y]=ode45(@rhs1,tspan,y0);
plot(t,y,'DisplayName',['y_0 = ',num2str(y0)],'linewidth',2)
hold on
end
title('directional field with seasonal effects: b = 0.005')
legend()
b = 0.05
figure(2)
yprime=@(t,y) 0.75.*y -0.05.*y.^2-1.5*y.^3./(1.25+y.^3).*abs(sin(pi.*t/365));
% Define the range of y and t values
y_min = 0;
y_max = 30;
t_min = 0;
t_max = 365*10;
% Define the number of grid points
num_points = 100;
% Generate grid points
yval = linspace(y_min, y_max, num_points);
tval = linspace(t_min, t_max, num_points);
hold on
dirfield(yprime,tval,yval)
for i=1:1:length(y0s)
y0=y0s(i);
[t y]=ode45(@rhs2,tspan,y0);
plot(t,y,'DisplayName',['y_0 = ',num2str(y0)],'linewidth',2)
hold on
end
title('directional field with seasonal effects: b = 0.05')
legend()
figure(3)
yprime=@(t,y) 0.75.*y -0.1.*y.^2-1.5*y.^3./(1.25+y.^3).*abs(sin(pi.*t/365));
% Define the range of y and t values
y_min = 0;
y_max = 10;
t_min = 0;
t_max = 365*10;
% Define the number of grid points
num_points = 100;
% Generate grid points
yval = linspace(y_min, y_max, num_points);
tval = linspace(t_min, t_max, num_points);
hold on
dirfield(yprime,tval,yval)
for i=1:1:length(y0s)
y0=y0s(i);
[t y]=ode45(@rhs3,tspan,y0);
plot(t,y,'DisplayName',['y_0 = ',num2str(y0)],'linewidth',2)
hold on
end
title('directional field with seasonal effects: b = 0.1')
legend()
table(t,y)
function [yp] = rhs1(t,y) % do not change this line
% *** Your entry 6 *** (start entering your system here)
a= 0.75;
b = 0.005;
p = 1.5;
q = 1.25;
yp=a*y -b*y.^2 -p*y.^3./(q+y.^3).*abs(sin(pi.*t/365));
end
function [yp] = rhs2(t,y) % do not change this line
% *** Your entry 6 *** (start entering your system here)
a= 0.75;
b = 0.05;
p = 1.5;
q = 1.25;
yp=a*y -b*y.^2 -p*y.^3./(q+y.^3).*abs(sin(pi.*t/365));
end
function [yp] = rhs3(t,y) % do not change this line
% *** Your entry 6 *** (start entering your system here)
a= 0.75;
b = 0.1;
p = 1.5;
q = 1.25;
yp=a*y -b*y.^2 -p*y.^3./(q+y.^3).*abs(sin(pi.*t/365));
end
Published with MATLAB® R2020a
