ODEs
From CDS 130
Contents 
1. Objectives
 To introduce the concept of an ordinary differential equation
2. Motivation
 Most science models have a mathematical representation that has the form of an ordinary differential equation.
3. Priming questions
 You won $100 and want to invest it. Which bank offers the better deal?
 Bank A: 12% interest per year compounded yearly
 Bank B: 1% interest per month compounded monthly
 Bank C: (1/365.2425)% interest compounded daily
"Compounded yearly/monthly/daily" means that at the end of one year/month/day the interest is added to your balance.
So, after one
 year you have $112 dollars
 month you have $101 dollars
 day you have 100+(1/365.2425) dollars
4. Notes
4.1. Introduction
In computational science, we most often want to simulate systems that continuously change. Previously we modeled population change and bank balance growth with a discrete equation. The mathematical model of this had the form
P(next instant)P(this instant) = a*P(this instant)
And the time between instances was some fixed time interval, say one month or one year.
4.2. Example
Population grows at a rate of 20% per year
P(next year) = P(this year) + 0.2*P(this year)
Will you get the same answer if you divide the growth rate by 12 and change the equation to go monthtomonth?
P(next month) = P(this month) + (0.2/12)*P(this month)
this can be rewritten as
P(next month)  P(this month) = (time between iterations)*(0.2)*P(this month)
where
(time between iterations) = (1/12)
4.3. Observation cont.
Will you get the same answer if you divide the growth rate by 12 and change the equation to go monthtomonth?
Answer: No. But, but the answers are close. We are going to use this fact to reduce the amount of computation that we need to do.
If we are told that the correct problem to solve is the one on the right, we can get an approximate answer that requires fewer steps by doing the problem on the left.
The difference between the answers is 1.2001.2194 = 0.006. If we are only concerned about getting the answer correct to two decimals, then it makes sense to solve the problem on the left if you want to do less computation.
The number of iterations is one and the percent increase per iteration is 20%. clear; P(1) = 1; for i = [1:1] P(i+1) = P(i) + (0.2/1.0)*P(i); end P(2) % Will display 1.2000 Longhand version of the above. clear; P(1) = 1; i = 1; P(i+1) = P(i) + 0.2*P(i); P(2) % Will display 1.2000 
The number of iterations executed is increased by a factor of 12 and the percentage increase in clear; P(1) = 1; for i = [1:12] P(i+1) = P(i) + (0.2/12)*P(i); end P(13) % Will display 1.2194 Longhand version of the above. clear; i = 1; P(i+1) = P(i) + (0.2/12)*P(i); i = 2; P(i+1) = P(i) + (0.2/12)*P(i); i = 3; P(i+1) = P(i) + (0.2/12)*P(i); i = 4; P(i+1) = P(i) + (0.2/12)*P(i); i = 5; P(i+1) = P(i) + (0.2/12)*P(i); i = 6; P(i+1) = P(i) + (0.2/12)*P(i); i = 7; P(i+1) = P(i) + (0.2/12)*P(i); i = 8; P(i+1) = P(i) + (0.2/12)*P(i); i = 9; P(i+1) = P(i) + (0.2/12)*P(i); i = 10; P(i+1) = P(i) + (0.2/12)*P(i); i = 11; P(i+1) = P(i) + (0.2/12)*P(i); i = 12; P(i+1) = P(i) + (0.2/12)*P(i); P(13) % Will display 1.2194 
4.4. Observation cont.
Consider what happens as we increase the number of iterations and decrease the multiplication factor by the same factor.
When the number of iteration is "large enough", you'll start getting the same final answer to several decimal places.
Importance: Most mathematical models imply that the number of iterations must be infinite, which would require and infinite amount of time to get an answer. In this example, we see that if we choose a large number, we'll get an answer that is expected to be close to the answer that we would get if we performed an infinite number of iterations.
clear; P(1) = 1; for i = [1:1] P(i+1) = P(i) + (0.2/1.0)*P(i); end P(2) % Will display 1.2000 clear; P(1) = 1; for i = [1:12] P(i+1) = P(i) + (0.2/12)*P(i); end P(13) % Will display 1.2194 clear; P(1) = 1; for i = [1:24] P(i+1) = P(i) + (0.2/24)*P(i); end P(25) % Will display 1.2204 
clear; P(1) = 1; for i = [1:100] P(i+1) = P(i) + (0.2/100)*P(i); end P(101) % Will display 1.2212 clear; P(1) = 1; for i = [1:1000] P(i+1) = P(i) + (0.2/1000)*P(i); end P(1001) % Will display 1.2214 clear; P(1) = 1; for i = [1:10000] P(i+1) = P(i) + (0.2/10000)*P(i); end P(10001) % Will display 1.2214 
4.5. Summary
In computational science, we most often want to simulate systems that have the computational model of (where The main problem is that if 
for i = 1:HUGE P(i+1) = P(i) + (a/HUGE)*P(i); end 
4.6. Definition
This type of program arises from the translation of a mathematical model that has the form of an ordinary differential equations (ODE) to a computational model. ODEs are used to mathematically model many natural systems.
for i = 1:HUGE P(i+1) = P(i) + (1/HUGE)*a*P(i); end
For example, the science model for a system that experiences continuous growth at a rate that is proportional to the population is an ordinary differential equation which has the form
 or
and Δt is a very small number (and takes the place of 1/HUGE
, which is a very small number). In this course, you don't need to worry about the mathematical notation  you should know that many natural systems follow an ordinary differential equation, and that you actually know how to solve ordinary differential equations on a computer!
4.7. Definition cont.
can be written as
P(next instant)P(this instant) = (time between this instant and the next)*a*P(this instant)
because by definition
and
The differential in differential equation comes from the fact that it is a model of a difference in population from one instant to the next for a difference in time that is very small.
4.8. Example
An example where an ODE appears is in the mathematical model of how an object cools. The change in temperature, T
, of an object as it cools has the mathematical model of
where T_{a} is the temperature of the room where the object is stored. This can also be written as
T(next instant)T(this instant) = (time between this instant and the next)*q*(T(this instant)Ta)
where the time between the next instant and this instant, Δt is nearly zero.
4.9. Summary
Any time you are asked to solve a problem that has this form:
P(next instant)P(this instant) = (time between this instant and the next)*a*P(this instant)
you should know that
 You are solving an ordinary differential equation.
 Even if the time between instants is very small (and so you need to do a large number of iterations to compute a final answer), you can get an approximate answer by choosing a time between instants that is "small enough" and a number of iterations that is "large enough".
5. Questions
5.1. Population Eqn.
Consider the equations
P(next year) = P(this year) + 0.01*P(this year)
P(next month) = P(this month) + (1/12)*0.01*P(this month)
 Write a program that uses the first equation to compute the answer after 10 years. What is P(11)?
 Write a program that uses the second equation to compute the answer after 10 years. What does this equation predict the population will be after 10 years have elapsed?
 Do 1. and 2. using a spreadsheet and write down your answers.
5.2. Population Eqn.
Consider the following model of population
Every year, population increases a value of 10% of the population in the previous year. However, if the predicted population is over 100, a disease outbreak instantly kills 80% of this predicted population value. For example, if the predicted population is 110, then the next year the population is 0.2*110.
 Use Excel to plot population as a function of time for 40 years. Assume that the initial population is 20.
 Try to make a change to the default setting for the plot that makes it clearer or easier to read (in your opinion; there is no "wrong" answer).
5.3. Heat Eqn.
The equation that describes how the temperature of an object, T
changes with time (Newton's Law of Cooling) is
 Use a spreadsheet to determine how the temperature of an object that starts out at 100 degrees C cools when the ambient temperature,
Ta
is 0. Assume that q = 0.1 and Δt = 0.1. At approximately what time does the temperature equal onehalf of its initial temperature?  Use a spreadsheet to determine if or how the time it takes for the temperature to fall by onehalf depends on the initial temperature of the object. Describe in words you how you used the spreadsheet to answer this question and your conclusion.
5.4. Code Example
Consider the following science model: "The starting balance in a bank account is $1,000 and it grows by 5% per year."
The following code was written to simulate the monthtomonth growth of this bank account for 36 months, and, to print out the balance at the end of months 26, 27 and 28:
B(1) = 1000; for i = [2:35] B(i) = B(i1) + 0.05*B(i1); end B(26) B(27) B(28)
When this code is run, we find that B(26) = $3,386.40, B(27) = $3,555.70 and B(28) = $3,733.50, all of which are not consistent with the science model. Why?
6. Activities
6.1. Interest
Suppose that you deposit $1000 in your account end of each year and that you started with $100 at the start of the first year.
Would you rather have a bank that offered 5% interest per year like this:
B(next year) = B(this year) + 0.05*B(this year) + 1000
or a bank that offered (5/12) percent interest per month like this:
B(next month) = B(this month) + (0.05/12)*B(this month) + 1000/12
Answer the following without using a computer program (work it out using pencil, paper, and a calculator). Turn in your work on a sheet of paper at the start of class (or post an image of your work take with a digital camera or scanner).
 For the first equation, write down the balance at the start of year one and year two.
 For the second equation, write down the balance at the start of the the start of each month up to an including month 25.
Write a computer program the does the above calculations. For both questions, a for
loop must be used.
Answer 

year 1 B = 100 year 2 B = 100 + 0.05*100 + 1000 = 1105 year 3 (note problem only asked for year 2, but year 3 included here for comparison with other answer) B = 1105 + 0.05*1105 + 1000 = 2160.25 B(1) = 100 for i = 1:2 B(i+1) = B(i) + (0.05/12)*B(i) + 1000/12 end B(3) month 1 B = 100 month 2 B = 100 + (0.05/12)*100 + 1000/12 = 183.3375 (Note  it is not correct to take 183.3375 and multiply it by 12! You must work out every iteration.) month 3 B = 183.3375 + (0.05/12)*183.3375 + 1000/12 = 267.4347 ... month 25 B = 2117.1659 + 2117.1659*(2117.16598644552) + 1000/12 = 2209.32 B(1) = 100 for i = 1:24 B(i+1) = B(i) + (0.05/12)*B(i) + 1000/12 end B(25) Note  the for loop approach should give the same answer as your hand calculation. Also note that the start of year 3 is month 25 which you computed. This is larger than the balance in all of year 3 using the first approach ($2160.25). Therefore, the second equation will result in more money in the account. 
6.2. Population Eqn.
Consider the equations
P(next year) = P(this year) + 0.01*P(this year)
P(next month) = P(this month) + (1/12)*0.01*P(this month)
Assume the initial population is 100.
 Write a program that uses the first equation to compute the population in year 10 and prints out its value.
 Write a program that uses the second equation to compute the population in the first month of year 10 and prints out its value.
 Write a program that uses the second equation to compute the population in the last month of year 10 and prints out its value.
When your solutions for each of the above are copied onto the command line, only the requested value must be printed. As an example, the following program prints out only the last value of an array created using a for
loop:
for i = [1:5] A(i) = i; end A(5)
Answer 

Note  only need to go till i = 9 to compute value in year 10. if you changed for i = [1:9] to for i = [1:10] and left everything else the same, you would get the correct answer, but you would have done an extra calculation that was not needed. clear; P(1) = 100; for i = [1:9] P(i+1) = P(i) + 0.01*P(i); end P(10) = 109.36 % 9 years x 12 months = 108. clear; P(1) = 100; for i = [1:108] P(i+1) = P(i) + (0.01/12)*P(i); end % Month 109 is first month of year 10 P(109) = 109.41 % 108 + 11 = 119 clear; P(1) = 100; for i = [1:118] P(i+1) = P(i) + (0.01/12)*P(i); end % Month 119 is last month of year 10 P(119) = 110.33 
6.3. Population Eqn.
Consider the following model of population
Every year, population increases a value of 10% of the population in the previous year. However, if the predicted population is over 100, a disease outbreak instantly kills 80% of this predicted population value. For example, if the predicted population is 110, then the next year the population is 0.2*110.
 Write a program that computes the population over a 40year time span. Assume that the initial population is 20. (Hint: you can do this using an
if
statement.)
Answer 

% Change 2400 to 39 to compute population over 40year time span. clear; P(1) = 20; for i = [1:2400] P(i+1) = P(i) + 0.1*P(i); if (P(i+1) > 100) P(i+1) = 0.2*P(i+1); end end plot(P) xlabel('Year') ylabel('Population')
