1. Finding Roots Algorithms
- Question 1.5: I chose to write the program in C. The code can be found here.
I investigated the behaviour of the Secant method with changes in the initial guesses of the root, and the variety of cases that were tested are described in the following table:
It can be noticed, that as long as one of the guesses is close to the root, there is pretty fast convergence, and that the closer the guess to the root, there is a faster convergence.
Also, it can be noticed, that the value of x1 is slightly more significant than the value of x0. But more importantly, since this is a simple quadratic function,
no matter how large the initial guess is, the algorithm converges to the correct answer. I even tried running it with x0=1000001 and x1=1000000 and received the answer within 32 iterations.
- Question 1.6:
This program is also written in C. The code can be found here.
When I ran the code, I saw that for initial guesses under 1 or very close to 1 the algorithm converges, but even for an initial guess of x=2 it diverges.
For x>>1 the method does not converge. This can be easily understood from the tanh graph - for x>>1 tanh(x) is approximately constant, meaning f'(x)=0 and therefore the
expression f/f' goes to infinity.
Now we wish to derive the critical value of the initial guess above which convergence will not occur.
This will happen if xi+1 is farther away from zero than xi, meaning: |xi+1|>|xi|.
Let's define g(x)=f(x)/f'(x), then we shall demand |xi-g(x)|>|xi|. The solution is g(x)>2x. Now, if x>x0, this directly implies that x2<-x1 and x1 < -x2 < x3 etc. The NR algorithm diverges.
In our case f(x)=tanh(x) and f'(x)=sec2(x), so the critical value of the initial guess is received by tanh(x)/sec2(x)=2x and the solution is x0=1.08866.
And indeed, for x0=1.08866 the code produces an answer of divergence but for x0=1.08865 we get the correct zero in 13 iterations.
I tried solving with the secant method as well (in the same code) and found that it converges for slightly higher initial guesses, but also has a limit.
For an initial guess of zero, the NR method converges in a single iteration. For the secant method, if x1 is zero it converges in a single iteration, not matter how large x0 is, but if x0 is zero and x1 is not it will take more than one iteration, depending on the x1 guess.
2. Ordinary Differential Equations
- Question 2.8:
I changed the potential to a central one: V(x,y)=2*sqrt(x2+y2)5. Following is the result for E=.1 and initial conditions y0=50, py0=2.
Please excuse me for not changing the titles of the graphs, I didn't want to mess with the code too much.
I tested this for many different initial conditions and found that if y0 is smaller the radius of the cross section of the toroid is bigger (infinity fr y0=0). For example, following is the same potential for y0=20 and the same py0=2.
Note that the radius is larger than before.
But, for the momentum it's the opposite - the larger py0 the larger the radius. For example, following is the same potential for pyy0=20 and the same y0=50.
In conclusion, the shape of the trajectory in the phase space does not change, but the initial conditions affect only the geometry a little. We recieve the toroidal shape, just as expected from Fig. 2.2 in the text book.
- Question 2.10:
For E=0.025 with time stamp 0.1 and initial conditions y0=2, py0=5, I recieved:
For E=0.05 and the same initial conditions, I recieved:
For E=0.075 and the same initial conditions, I recieved:
For E=0.1 and the same initial conditions, I recieved:
For E=0.15 and the same initial conditions, I recieved:
In all graphs elliptic fixed points can be found in the centers of the slliptic surfaces. chaotic regions can hardly be found. Therefore, I tried changing the initial xonditions.
For E=0.15 and larger y0, I recieved:
This is still not chaotic.
For E=0.15 and larger py0, I recieved:
Finally we recieve a chaotic system.
3. Using the Runge-Kutta Method to solve the Hodgkin-Huxley Model
I choose to describe a problem discussed in an undergraduate course I have taken in the EE department - Introduction to Biological Signals and Systems.
In physiology, an action potential is a short-lasting event in which the electrical membrane potential of a cell rapidly rises and falls, following a consistent trajectory.
The shape of a typical action potential is depicted in the following figure. The membrane potential remains near a baseline level until at some point in time it abruptly spikes upward and then rapidly falls.
The Hodgkin–Huxley model, or conductance-based model, is a mathematical model that describes how action potentials in neurons are initiated and propagated, through a system of four first-order non-linear differential equations:
where V is the membrane potential, I is the externally applied current. VNa
are the ions' resting potential. The current of each ions type is proportional to the potential difference between the membrane potential and the rest potential. The parameters gi
and the rest potentials are dependent on the neuron type and are found experimentally. Following is an example set of parameters:
We wish to solve the above set of equations for V(t), given a set of initial conditions and a given externally applied current.
I would chose to solve this problem using the fourth order Runge-Kutta method. This method is chosen because it is more accurate than the Euler method and since the equations are not linear, the implicit methods can be difficult.
The steps for solving would be the following:
- Substitute the initial V(1) to calculate α(V(1)),β(V(1)) and substitute to get m(1),n(1),h(1).
- For each time step i: we shall define ni=dn/dt(i) and the same for mi and hi and calculate
- k1=dt*HH(i,V(i),n(i),m(i),h(i)) where HH is the differential equation for dV/dt.
- Create the next step for each variable:
I would use MATLAB since it is the easiest and the language I am used to work with, but it can be easily implemented in any language.
For a given set of initial conditions and not too long a time vector, this should not be too long of a calculation, and I would plot V(t) at the end.
But, if we wish to examine the response to changing external currents for a long period of time, I might choose to plot the progress of V(t) at every iteration.
Accuracy and Validation
The accuracy is the order of dt4
, so the smaller dt the better. But, as we've learned, this can be a problem if dt is smaller than the precision of the arithmetic in the computer.
Another problem might occur if the initial conditions and external current are such that will create chaos. Reduction of the above equations to two dimensions under certain assumptions can be examined in state-space representations and there isn't always a steady state. In real life this cannot happen (since n and h bring V down) but if the given conditions are not realistic one might receive chaos.
This should not happen.