Difference between revisions of "Physics 212, 2019: Lecture 6"

From Ilya Nemenman: Theoretical Biophysics @ Emory
Jump to: navigation, search
(Created page with "{{PHYS212-2019}} ==An example of a model: Solving for an equilibrium position, harder version== ;Problem formulation: A charged particle free to move on a rod is put in the...")
 
 
(6 intermediate revisions by the same user not shown)
Line 1: Line 1:
 
{{PHYS212-2019}}
 
{{PHYS212-2019}}
  
−
==An example of a model: Solving for an equilibrium position, harder version==
+
==Introduction to verification of models==
 +
The most important thing to remember here is that a computer does exactly what you ask it to do. And we are pretty bad at following our own instruction in an unbiased manner. You will be reading the code and thinking that it does what you want it to do, but you will be wrong. Your biases and expectations will have the best of you. Thus it is almost useless to verify the code line by line for correctness (unless we are dealing with a simple syntactic error). Instead one needs to run the code many times, to look at what it outputs, and to ask if the output makes sense.
  
−
;Problem formulation:
+
This is not very precise -- and, indeed, verifying a computational model is often as much an art as it is science. But some simple rules of thumb apply. First, ''do not verify your model all at once!'' Instead, just like we build models hierarchically, we need also to verify them hierarchically. For this, you verify and seal every leaf on the tree that is your model. Then you verify blocks that consist of many leafs, then the blocks of blocks, and so on. Do not attempt to write (and certainly not to verify) a larger piece of a model when its constituent pieces are not yet verified to your satisfaction.
−
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: The analysis above suggests that the entire model consists of solving for 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>. This is the first of many computational algorithms that we will introduce in this class, and it allows to find roots (zero crossings) of more-or-less arbitrary functions, not just of polynomials.
 
  
−
:Let's suppose that we have a guess <math>x_0</math> for a root of a function <math>F(x)</math>. 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 one can also choose to interrupt the solution choice when the function itself is approximately zero, and both approaches are valid. They are, in fact, equivalent for typical functions, but may give different results for various special cases -- and the discussion of this in "Numerical recipes" is illuminating.  
+
Second, when verifying an piece of the code, every free parameter, every variable in the code must be verified. For this, you need to set every one of the parameters (and some parameter combinations) in turn to values that make the problem easier. These could be value that make your problem analytically solvable, or reduce your problem to something that you already solved before. And then you must ensure that the output of your code with such ''special cases'' parameter values is exactly what you expect it to be from the previous solutions. There's something called the thirteenth strike rule -- if a clock strikes 13 times, then all of them (not just the thirteenth) are wrong. If you see that the problem does something unreasonable for these special cases, it's more likely than not that it is also wrong in general. If you found a mistake -- fix it.
  
−
:A few notes about the algorithm are in order. First, the algorithm will only find a root if it is initialized with the initial guess in the root's 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, and then run the algorithm for each of the guesses. Even then the solution may not exist -- there is no guarantee that the the iterations converge; they can enter cycles or diverge instead. Most of such problems can be tracked to the fact that the iteration step involves dividing by <math>G</math>, derivative of the function <math>F(x)</math>. When this derivative becomes small, changes to <math>x</math> will be large, and all sorts of problems will emerge. We will explore this in the ''your work'' section.
+
Third, spend time on verification and try all sorts of initializations, even if you do not expect that the initialization simplifies the problem, allowing you to compare to the prior knowledge. Instead, what is likely to happen is that, for some of the initial conditions, your model will produce obviously wrong results. For example, it may crash, or may produce NaN outputs, or will enter an infinite loop, or produce negative numbers for physical quantities that can only be positive, or something similar. And once you see a problem, do not just step over it -- solve it! Figure out what's going on, and fix things by either introducing more exception checks, or by fixing the logic of the solution. The more things you try, the likely you are to detect more errors. But also remember -- small changes in parameters are almost useless. Go for large changes, explore the edges of the parameter spaces. You are more likely to identify problems this way.
  
−
:'''Your work:''' with the provided skeletal 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? If so, do the modification. Can you understand why some of your attempts converged to a solution, and some did not? '''Don't forget to submit your code.'''
+
The end of the [[Physics 212, 2019: Lecture 5|Lecture 5]] has an example of how one would verify a computational model, namely the model of finding the equilibrium position of a charge surrounded by three other charges. Much of the course this semester is about how to verify models -- and as you practice this, you will become better at this task.
−
 
−
:Once the Newton solver for roots has been implemented, the rest of the model design is simple. We start with initializing variables, verifying various possible exceptions, finding the root of the force equal to zero equation, and then outputting the root. A script for this is available for download,
 
−
as a [[media:Newton2019.txt | Jupyter notebook]].
 
−
 
−
;Model verification: To verify the model, we need to check that it produces results that we know to be correct for various special cases. As a general rule, every term, every function in the model must be verified. For this, you set all of the parameters in turn to values that make the problem easy, and compare the output to what you should have expected. In this particular case, the following verifications would be reasonable:
 
−
:*Set the initial condition outside of the leftmost or the rightmost charge. Then there should be no equilibrium, and the charge will move to infinity -- which will mean that the Newton solver won't converge. Do we observe this?
 
−
:*If we have three fixed charges, but one of them is zero, then the problem is equivalent to having two fixed charges only. And that problem we already solved previously. So we can set, in turn, the values of each of the three fixed charges to zero, and compare that the output of the program is the same as that of the two-charges quadratic solver program. An even simpler first step could be setting one of the charges to zero, and the two other equal to each other. Then the movable charge should settle exactly in the middle between the two non-zero charges. Does this?
 
  
−
: Notice that with these checks, we verify the effects of the initial condition and the value and the position of the fixed charges to be what we expect them to be. Thus the verification is more or less complete. In addition, none of the verification steps is useless -- if you only checked the initial condition, then the effect of the fixed charges would not be explored, and so on.
+
==Jupyter notebooks and writing reports==
−
+
Appendix B in the textbook has a nice introduction of what Jupyter notebooks are. They are basically lab notebooks for computational experimentation and are a great way of writing lab reports. They allow for easy embedding and execution of Python code, of equations, of external pictures, and of Python graphics output.  
−
;Discussion: Here we implemented a solution for the equilibrium position of the charge in the background of three fixed charges.  The program, based on the Newton-Raphson algorithm, found the solution. Exceptions for position of the charges and the charge values were encoded. The program output has been verified against solution of a simpler problem with only two fixed charges, and the results (hopefully) matched. We assumed here that only the equilibrium position mattered. In principle, we could have written the problem as a dynamical model, and used the Newton's laws to find the position of the charge as a function of time. The value of the position at large time would be the equilibrium position, provides there is some friction in the problem that would force the charge to stop eventually. Another simplification we assumed was that the charges are positioned in just one dimension (on a rod). For a 2-d or a 3-d configuration of the charges, we would need to extend the Newton-Raphson algorithm to be useable for multi-dimensional root finding, which goes slightly outside of the scope of this class (would require knowing multivariate calculus).
 
  
−
==New Python Concepts==
+
I provided for you a few sample Jupyter notebooks, including
−
*Loops of different types (for and while).
+
*[[media:Newton2019.txt | Jupyter notebook for solving the equilibrium position of a charge problem]].
−
*We started talking about scopes: what happens to the variable ''i'' after the ''for'' loops finishes executing? 
+
*[[media:SampleReport2019.txt | Sample project report in a Jupyter notebook]] for the problem referenced above.
−
*Variables vs objects, yet again.
 
−
*Vector math is always faster than element-by-element math -- we showed this in a simple script in class.
 
  
−
==Old Math Concepts==
+
There are many online resources about Jupyter notebooks -- and I suggest you read some of them (see below for one link). These notebooks may be a very useful tool for you in the future. For example, in my own research group, students write weekly progress reports and presentations for group meetings using Jupyter notebooks.
−
*Make sure you review Taylor series. These will be useful in the rest of the course.
+
*[https://github.com/jupyter/jupyter/wiki/A-gallery-of-interesting-Jupyter-Notebooks#introductory-tutorials Introductory tutorials for Jupyter notebooks]
 +
 
 +
==Quiz 1==
 +
At the end of the lecture, Quiz 1 will open up on Canvas, and you must finish it and submit by the end of the class.

Latest revision as of 15:16, 3 February 2019

Emory Logo

Back to the main Teaching page.

Back to Physics 212, 2019: Computational Modeling.

Introduction to verification of models

The most important thing to remember here is that a computer does exactly what you ask it to do. And we are pretty bad at following our own instruction in an unbiased manner. You will be reading the code and thinking that it does what you want it to do, but you will be wrong. Your biases and expectations will have the best of you. Thus it is almost useless to verify the code line by line for correctness (unless we are dealing with a simple syntactic error). Instead one needs to run the code many times, to look at what it outputs, and to ask if the output makes sense.

This is not very precise -- and, indeed, verifying a computational model is often as much an art as it is science. But some simple rules of thumb apply. First, do not verify your model all at once! Instead, just like we build models hierarchically, we need also to verify them hierarchically. For this, you verify and seal every leaf on the tree that is your model. Then you verify blocks that consist of many leafs, then the blocks of blocks, and so on. Do not attempt to write (and certainly not to verify) a larger piece of a model when its constituent pieces are not yet verified to your satisfaction.

Second, when verifying an piece of the code, every free parameter, every variable in the code must be verified. For this, you need to set every one of the parameters (and some parameter combinations) in turn to values that make the problem easier. These could be value that make your problem analytically solvable, or reduce your problem to something that you already solved before. And then you must ensure that the output of your code with such special cases parameter values is exactly what you expect it to be from the previous solutions. There's something called the thirteenth strike rule -- if a clock strikes 13 times, then all of them (not just the thirteenth) are wrong. If you see that the problem does something unreasonable for these special cases, it's more likely than not that it is also wrong in general. If you found a mistake -- fix it.

Third, spend time on verification and try all sorts of initializations, even if you do not expect that the initialization simplifies the problem, allowing you to compare to the prior knowledge. Instead, what is likely to happen is that, for some of the initial conditions, your model will produce obviously wrong results. For example, it may crash, or may produce NaN outputs, or will enter an infinite loop, or produce negative numbers for physical quantities that can only be positive, or something similar. And once you see a problem, do not just step over it -- solve it! Figure out what's going on, and fix things by either introducing more exception checks, or by fixing the logic of the solution. The more things you try, the likely you are to detect more errors. But also remember -- small changes in parameters are almost useless. Go for large changes, explore the edges of the parameter spaces. You are more likely to identify problems this way.

The end of the Lecture 5 has an example of how one would verify a computational model, namely the model of finding the equilibrium position of a charge surrounded by three other charges. Much of the course this semester is about how to verify models -- and as you practice this, you will become better at this task.

Jupyter notebooks and writing reports

Appendix B in the textbook has a nice introduction of what Jupyter notebooks are. They are basically lab notebooks for computational experimentation and are a great way of writing lab reports. They allow for easy embedding and execution of Python code, of equations, of external pictures, and of Python graphics output.

I provided for you a few sample Jupyter notebooks, including

There are many online resources about Jupyter notebooks -- and I suggest you read some of them (see below for one link). These notebooks may be a very useful tool for you in the future. For example, in my own research group, students write weekly progress reports and presentations for group meetings using Jupyter notebooks.

Quiz 1

At the end of the lecture, Quiz 1 will open up on Canvas, and you must finish it and submit by the end of the class.

Retrieved from "https://nemenmanlab.org/~ilya/index.php?title=Physics_212,_2019:_Lecture_6&oldid=3177"