Physics 212, 2019: Lecture 4 - Revision history

Ilya at 18:55, 3 February 2019

2019-02-03T18:55:59Z

Ilya: Created page with "{{PHYS212-2019}} In this lecture, we discuss the basics of the modeling process, focusing on simple linear (exponential) bacterial growth as an example. ==What is a model?==..."

2019-02-03T18:11:22Z

Created page with "{{PHYS212-2019}} In this lecture, we discuss the basics of the modeling process, focusing on simple linear (exponential) bacterial growth as an example. ==What is a model?==..."

New page

{{PHYS212-2019}}

In this lecture, we discuss the basics of the modeling process, focusing on simple linear (exponential) bacterial growth as an example.

==What is a model?==
The first question we need to address is: What is a model?

;Model: A model in science is a simplified representation of a phenomenon, a system, or a process being studied.

Typically we build models because their simplified structure makes it easier for us to comprehend them, to analyze them, and to make predictions about how they will respond in yet un-tested scenarios compared to doing the same for the true system. In many respects, everything that you have learned about the scientific method in high school is wrong, or, at best, misleading: science is not about hypothesis testing. Science is about building, verifying, and improving models of Nature. (There is a great short recent article on this topic, which I encourage you to read: [https://www.wired.com/2017/01/dont-dare-try-teach-science-without-building-models/]). A globe is a model of the Earth. A double helix is a model of DNA. Newton's Second Law is a model of how material bodies affect each other. And a mouse, a fly, or an yeast are models of various aspects of human biology.

What is a good model? There's no unique answer to this question. Is Newton's Second Law a good model of interactions among bodies? It's an absurdly good model for the world we experience daily, which exists on scales of meters, kilograms, and seconds. But as soon as velocities become large (comparable to the speed of light), or masses become atomic or smaller, then the Second Law ceases being a good model, and relativistic equations of motion, or Schroedinger's equation of quantum mechanics become better. Similarly, a mouse is a good model of a human when we talk about basic cellular processes, but few would argue that it is a good model of higher level human cognition (indeed, unlike you, no mouse outside of the ''Hitchhiker's Guide to The Galaxy'' Universe will ever be able to read and comprehend this text). The upshot of this is that the quality of the model is determined by the question this model was designed to answer. The same model can be good in one context, and bad in another. The quality of the model is thus determined by comparing its predictions to those of experiments on a ''real system'' within the specific context of questions that this model is supposed to answer.

There are different kinds of models: physical models, material models, animal models, conceptual models, and many others. For example, a double helix is a model of DNA. A mouse is a model of some physiological processes in a human. However, in this class, we will talk about '''mathematical'' and ''computational'' models only, but all of the considerations above (and many from below) apply to many other kinds of models as well.

==Steps in the computational model-building process==
;Analysis of a problem
:Here, by reading assignments, literature, or talking to your colleagues/users, you get answers to the following questions. What is the question being asked? Why is this an interesting question? What is known about the problem and the answer? What form is the answer expected to be in? You often need to slightly rephrase the question being asked as a result of this analysis, making it more precise and focused. You also figure out which kind of information you need in order to solve the problem, and what is missing.

;Model development/formulation
:This step involves translating the problem into the actual mathematical model. It consists of the following sub-steps.
:*Gathering relevant data - Having analyzed the problem (above), we know which data are needed to be able to formulate it. These data need to be gathered from the prior literature, from experiments, or from other sources.
:*Listing and substantiating assumptions - Models are always ''simplified representations'' of the real processes or systems, and, therefore, are by definition wrong. These simplified assumptions must be listed explicitly, so that we know when to expect the model to be drastically wrong (if the assumptions are violated), or approximately correct (if the assumptions are not violated).
:*Determining variables and their units - Here we list all of the variables that are either dynamical (changing) or constant, and specify their units. Unit specification is important. Recall the story of the [https://en.wikipedia.org/wiki/Mars_Climate_Orbiter#Cause_of_failure Mars Climate Orbiter], which disintegrated in the Martian atmosphere because of inconsistent specification of units of measurements.
:*Determining relation among variables - Some variables will be ''constant'', and these need to be specified. Other variables depend on the rest by means of algebraic relations, and we will call them ''dependent'' variables. Yet other variables dynamically change, so that their current state depends on their past state and the current state of other variables; we will call these ''dynamical'' or ''container'' variables by analogy with a physical container, amount of material in which depends on the instantaneous flow rate and the amount of material in the past. Specification of dependencies among variables is typically done in charts -- we will denote constants as triangles, dependent variables as circles, and containers by boxes. Relations between the variables are denoted by arrows.
:*Writing down equations or rules - Finally, having identified all dependencies, we need to put a mathematical law at each arrow, identifying ''precisely'' how variables depend on each other

;Model implementation
:For computational models, this step involves writing a program to solve the problem using your favorite computer language (Python for our course).
:*Writing the model in your computing language of choice - Following the discussion on algorithmic thinking, we first write down an algorithm for solving the problem, and write down implementation of the algorithm in the computer language that we use. Needless to say, here we also verify that our program actually works -- that is, executes on a computer.
:*Solving the model using the tools of the language - Different languages will have different pre-programmed capabilities, implemented by previous generations of software developers. It is important to realize which tools are available to us and use such pre-developed tools in our implementation.

;Model verification
:This is a crucial step in the modeling process, which is often not discussed explicitly in many textbooks. As anecdote, I point out that almost every faculty member will be able to share a story with you, which will go roughly as follows. A student spends weeks on coding a solution to a certain problem, and then s/he comes to the professor with the result. It takes the professor just one question, one run of the program to show that the solution is wrong. And the student then leaves frustrated that s/he had spent so much time with nothing to show for it, and feeling very much down about himself or herself because of how easy it seemed it was for the professor to solve the problem. In fact, the professor didn't solve the problem. S/he just found a mistake in the student's solution, which was rather easy to do because the solution had not been tested/verified.
:*We verify the correctness of the solution first by verifying our code line by line, then block by block. But, crucially, we also verify it by testing special cases. For every interesting parameter/variable/interaction in the problem, there are always special cases of values of the corresponding parameters where the problem becomes easy (or at least easier) to solve, sometimes even analytically. So one sets the parameter to the special value and verifies if the program outputs the simple solution that we know it should. If it doesn't -- the program is wrong. This needs to be repeated for every important parameter in the problem before one can conclude that the solution is probably verified and is probably correct. The more independent special cases are verified, the higher is the probability that there is not a mistake in the solution.

;Interpretation and reporting of the results
This part will change depending on the specifics of the problem you are solving. However, generally, it involves:
:*Making plots, tables, or other visualizations of the program output.
:*Responding to the main question, for which the program was designed.
:*Discussing if the found solution is what we expected, and why or why not.
:*Discussing and interpreting the physical meaning of the solution.
:*Discussing what would happen to the solution if some of the simplifying assumptions get relaxed.

Finally, your report of the problem solution should follow the same steps as the modeling process and should contain '''all the same sections'''.

==Types of models==
There are a lot of different computational models that we will be exposed to in the course of this class. And there are even more that we won't be. Some specific types for you to keep in mind are the following:

;Probabilistic (stochastic) vs. deterministic models
:A probabilistic model is the one whose solution involves an element of chance, so that, even if run with the same conditions, detailed solutions of the model might be different in different runs. In contrast, a deterministic model does not have an element of chance, and so the solution is always the same if run with the same initial conditions.

;Static vs. dynamic models
:Static models are such that the variables we study do not depend on time. In dynamic models, the variables of interest depend on time.

;Spatially extended vs. point (or well-mixed) models
:In spatially-extended models, the solution is given by variables (field, as we will call them later) that are different at every point in space. In well-mixed models, only a small set of variables, independent of the spatial coordinates, characterizes the solution

;Continuous time vs. discrete time models
:In continuous time dynamic models, time changes continuously (though on a computer, time is always measured in discrete chunks). In discrete time models, time changes in specific, well-separated steps.

;Continuous space vs. discrete space models
:Similarly to the above, a spatially extended is discrete if the space is specified as a lattice, and it's continuous if every point in space can be considered.

Which specific models to use depends on the type of problem you study, and the choice, the explanation of it, and the assumptions involved, should figure prominently in the model-building process.

==An example of a model: Solving for an equilibrium position, light version==
A charged particle free to move on a rod is put in the between two particles with an electric charge of the same sign affixed to the end of the rod. Where will the charge come to rest?

;Analysis: How should we model this? is the model dynamic or static? At the first sign, this is a dynamic model, since we are talking about ''motion'' of a charge. On the other hand, we are asked only about the final equilibrium position. Then maybe this is a static model after all -- we only need to know where the charge comes to rest. The model is definitely a continuous space model -- position of the particle is a continuous variable. The model is also deterministic -- there's not a stochastic element present. To make this a static model (so that there's a final rustic point for the charge), there has to be friction, so that the particle does not keep oscillating infinitely. Nothing is mentioned about this in the problem formulation -- but it seems like this is presumed.Does it even make sense that the particle will come to a rest? Same sign particles repel. And the repulsion is stronger when they come closer. So a particle between two other fixed one will be pushed away from both ends -- it seems reasonable that it will come to a steady point at some time if friction robs it of energy.

;Model building: The equilibrium for the system will be <math>\sum \vec{F}=0</math>. Since the particle is moving on a rod -- in one dimension -- this is equivalent to <math>\sum {F}=0</math>. There are two Coulomb forces acting on the particle from the two end charges. Let's say the charges are at positions <math>x_1</math> and <math>x_2</math>, and have charges <math>q_1</math> and <math>q_2</math>respectively. Let's suppose <math>x_1<x<x_2</math>, where <math>x</math> is the position of the moving charge. The forces acting on the moving charge, which we call <math>q</math> are then <math>\sum F= \frac{qq_1}{(x-p_1)^2}-\frac{qq_2}{(x-p_2)^2}</math>. Thus to solve the problem, we need to solve the following equation for <math>x</math>: <math>\frac{qq_1}{(x-p_1)^2}-\frac{qq_2}{(x-p_2)^2}=0</math>, or, cancelling <math>q</math>, <math>\frac{q_1}{(x-p_1)^2}-\frac{q_2}{(x-p_2)^2}=0</math>. Interestingly, the only variables we need for initialization are the three charge values, and the two end charge positions -- the initial position of the moving charge, and charge value itself seem to be not needed. Multiplying the equation by the two denominators (which, luckily, cannot be zero), we get: <math>{q_1}{(x-p_2)^2}-{q_2}{(x-p_1)^2}=0</math>. Simplifying the equation, we get: <math>(q_1-q_2) x^2 +2(-q_1p_2+q_2p_1)x+(q_1p_2^2-q_2p_1^2)=0 </math>. So it seems that solving this problem is equivalent to solving a simple quadratic equation.

;Model implementation: To build a model, we start with large blocks describing chunks of the solution of the problem. We progressively break them apart into smaller and smaller blocks, until we are able to write the code for each one of the blocks. At the same time, as we are implementing the blocks, we collect pieces that need to be initialized in the initialization block, and need to be reported in the termination block. Here we realize that, in the previous class hour, we wrote a simple program that solves the quadratic equation <math>ax^2+bx+c=0</math>. Thus implementing the model is equivalent to taking the previously written solution for the quadratic. Thus the implementation will consist of just three blocks: (1) initialization, where we assign <math>a=q_1-q_2</math>, <math>b=2(-q_1p_2+p_2x_1)</math>, and <math>q_1p_2^2+q2p_1^2</math>; (2) solving the quadratic equation and getting its two solutions; and (3) termination, where, of the two solutions, we then output the one that falls between <math>p_1</math> and <math>p_2</math>.

;Model verification: To verify that the solution is correct, we should try it in cases where the solution is easy to guess or to estimate otherwise, and compare. For example, we can try <math>q_1=q_2</math>, where it's clear that the solution should be <math>x=(p_1+p_2)/2</math>. We could also choose a few conditions with different <math>q_1</math>, and we would expect that growing <math>q_1</math> should push the equilibrium farther from the left end.

;Discussion: Having solved the problem, we would usually discuss our findings here, and reiterate various constraints and assumptions that entered the solution. However, this problem is too simple too write much here.

;Your work: Implement the code for solving this problem, using your previous solution of the quadratic equation. Think about various exceptions that can happen with the solution. '''Make sure to submit your code.'''

==An example of a model: Solving for an equilibrium position, harder version==
A charged particle free to move on a rod is put in the between the first and the second of three immovable particles with an electric charge of the same sign affixed to the rod at various points. Where will the charge come to rest?

;Analysis: The same analysis as in the simpler version applies here: it's a static, deterministic, continuous space model. We again need the friction for the solution to exist, and the moving particle with come to rest, and the solution will exist.

;Model building: As before, the equilibrium for the system will be given by <math>\sum {F}=0</math>. There are now three Coulomb forces acting on the particle from the fixed charges. Let's say the charges are at positions <math>p_1<x_2<p_3</math>, with charges <math>q_1</math>, <math>q_2</math>, and <math>q_3</math>, respectively. Let's suppose again that the position of the movable particle is <math>x_1<x<x_2</math>. Then the solution of the problem will be given by <math>\frac{q_1}{(x-p_1)^2}-\frac{q_2}{(x-p_2)^2}-\frac{q_3}{(x-p_3)^2}=0</math>. We notice again that the three charges and fixed charge positions are the only things that are needed for the solution. This is equivalent to <math>{q_1}{(x-p_2)^2}(x-x_3)^2-{q_2}{(x-p_1)^2}(x-p_3)^2 -{q_2}{(x-p_2)^2}(x-p_1)^2=0</math>. This is not a quartic equation in <math>x</math>, and solving this problem would be equivalent to finding the root of this quartic equation between <math>p_1</math> and <math>p_2</math>. This is going to be more interesting!

;Model implementation: We need to find roots of a fourth order polynomial. How do we do this? For this we will introduce the Newton-Raphson method for finding a root of an equation <math>F(x)=0</math>. Let's suppose that we have a guess <math>x_0</math> for a root. We evaluate the function, <math>F_0=F(x_0)</math>, and it turns out that the guess is not quite right, <math>F_0\neq 0</math>. We can then write the Taylor expansion for the function <math>F</math> for the values of <math>x</math> around <math>x_0</math>. Namely, <math>F(x)\approx F_0 + G_0(x-x_0)</math>, where <math>G_0</math> is the value of the derivative of <math>F</math> evaluated at <math>x_0</math>, <math>G_0=\left(dF/dx\right)_{x_0}</math>. But then we can get a better estimate of where the root is by solving the equation <math> F_0+G_0(x-x_0)=0</math>. This would give <math>x_1=x_0-F_0/G_0</math> for the next estimate. Iterating this, one will get <math>x_{i+1}=x_i -F_i/G_i</math>. This way we will approach the solution progressively better and better, as we iterate more. At some point we will need to stop the iteration. Usually, this is yet another parameter that one adds to the initialization in the experiment: the tolerance <math>\epsilon</math> for finding the solution. That is, one interrupts the recursion when <math>\left|x_{i+1}-x_i\right|<\epsilon</math>. Note that the algorithm will only find a root if it is initialized with the initial guess in its vicinity. So if one wants to find a root, one needs to have good guesses. And to find all roots, one needs to have initial guesses near all of the roots. Even then the solution may not exist -- there is no guarantee that the the iterations converge; they can enter cycles or diverge instead. '''Your work:''' with the provided code implementing the Newton-Raphson method, try to find '''all''' roots of the following functions <math>F(x)= {\rm atan} 3x</math> and <math>F(x)=\frac{1}{3+x^2}-\frac{1}{6}</math>. Do you need to modify the code to deal with convergence of the iterations for some conditions? Can you understand why some of your attempts converged to a solution, and some did not? '''Don't forget to submit your code.'''

''More to be added.''

;Model verification:
;Discussion:

@@ Line 1: / Line 1: @@
 {{PHYS212-2019}}
 ==An example of a model: Solving for an equilibrium position, light version==
@@ Line 80: / Line 16: @@
 ;Your work: Implement the code for solving this problem, using your previous solution of the quadratic equation. Think about various exceptions that can happen with the solution. '''Make sure to submit your code.'''
-==An example of a model: Solving for an equilibrium position, harder version==
+==New Python concepts==
-A charged particle free to move on a rod is put in the between the first and the second of three immovable particles with an electric charge of the same sign affixed to the rod at various points. Where will the charge come to rest?
+In this class, we introduced and discussed the following concepts:
+*Objects -- mutable (arrays and lists) vs. immutable (numbers)
-;Analysis: The same analysis as in the simpler version applies here: it's a static, deterministic, continuous space model. We again need the friction for the solution to exist, and the moving particle with come to rest, and the solution will exist.
+**Object attributes and object methods, using dir()
+**Overloading methods
-;Model building: As before, the equilibrium for the system will be given by <math>\sum {F}=0</math>. There are now three Coulomb forces acting on the particle from the fixed charges. Let's say the charges are at positions <math>p_1<x_2<p_3</math>, with charges <math>q_1</math>, <math>q_2</math>, and <math>q_3</math>, respectively. Let's suppose again that the position of the movable particle is <math>x_1<x<x_2</math>. Then the solution of the problem will be given by <math>\frac{q_1}{(x-p_1)^2}-\frac{q_2}{(x-p_2)^2}-\frac{q_3}{(x-p_3)^2}=0</math>. We notice again that the three charges and fixed charge positions are the only things that are needed for the solution. This is equivalent to <math>{q_1}{(x-p_2)^2}(x-x_3)^2-{q_2}{(x-p_1)^2}(x-p_3)^2 -{q_2}{(x-p_2)^2}(x-p_1)^2=0</math>. This is not a quartic equation in <math>x</math>, and solving this problem would be equivalent to finding the root of this quartic equation between <math>p_1</math> and <math>p_2</math>. This is going to be more interesting!
+*Distinction of variables and objects; we looked at a few examples of operations that create new objects, and a few that do not.
+*We also discussed handling exceptions -- such as when the discriminant in the quadratic equation is zero, or when roots are complex.
-;Model implementation: We need to find roots of a fourth order polynomial. How do we do this? For this we will introduce the Newton-Raphson method for finding a root of an equation <math>F(x)=0</math>. Let's suppose that we have a guess <math>x_0</math> for a root. We evaluate the function, <math>F_0=F(x_0)</math>, and it turns out that the guess is not quite right, <math>F_0\neq 0</math>. We can then write the Taylor expansion for the function <math>F</math> for the values of <math>x</math> around <math>x_0</math>. Namely, <math>F(x)\approx F_0 + G_0(x-x_0)</math>, where <math>G_0</math> is the value of the derivative of <math>F</math> evaluated at <math>x_0</math>, <math>G_0=\left(dF/dx\right)_{x_0}</math>. But then we can get a better estimate of where the root is by solving the equation <math> F_0+G_0(x-x_0)=0</math>. This would give <math>x_1=x_0-F_0/G_0</math> for the next estimate. Iterating this, one will get <math>x_{i+1}=x_i -F_i/G_i</math>. This way we will approach the solution progressively better and better, as we iterate more. At some point we will need to stop the iteration. Usually, this is yet another parameter that one adds to the initialization in the experiment: the tolerance <math>\epsilon</math> for finding the solution. That is, one interrupts the recursion when <math>\left|x_{i+1}-x_i\right|<\epsilon</math>. Note that the algorithm will only find a root if it is initialized with the initial guess in its vicinity. So if one wants to find a root, one needs to have good guesses. And to find all roots, one needs to have initial guesses near all of the roots. Even then the solution may not exist -- there is no guarantee that the the iterations converge; they can enter cycles or diverge instead. '''Your work:''' with the provided code implementing the Newton-Raphson method, try to find '''all''' roots of the following functions <math>F(x)= {\rm atan} 3x</math> and <math>F(x)=\frac{1}{3+x^2}-\frac{1}{6}</math>. Do you need to modify the code to deal with convergence of the iterations for some conditions? Can you understand why some of your attempts converged to a solution, and some did not? '''Don't forget to submit your code.'''
+*Code re-use. We re-used the quadratic equation solution code from the previous lecture here to find the solution for the charge position.