https://docs.google.com/a/vassar.edu/document/d/1nvyB5giUicIIiCJuKNKTwJagBlepP8J-IOQsMeMgB7Q/edit

Brownian Motion of a Drunkard

Random Unit Vector 2D Motion of a Drunkard

As shown by the standard error of the slope within slopes of the displacement squared, the Brownian Motion has a statistically significant larger displacement squared. This is before I alter the 2D code to account for more drunken chance variables. This will serve as a good baseline to see how the mean displacement changes after each variable is added or combinations of variables are tested. I chose to use brownian motion as the standard since it is purely random in every way (in my code). The direction of a step and length are always randomized. I am very pleased with how my codes came out and the visualizations they provide. I apologize for the poor quality of screenshots but they will not look this way in the final project. The top subplot is one representation of the random movement. The second is a displacement squared of one walk. The final in blue is the mean of all of the displacement squareds. The red is the best fit line with the statistics below.

This will serve as my jumping off point to compare Brownian motion of a drunkard (a gas) to a college student on the weekend with several variables and options besides just taking steps. I am curious to see how much I can alter the mean displacement squared with various variables and different weights. Unfortunately, I still need to collect more data outside of class to correctly weight the variables. I’m still catching up on other work from tech week.

After investigating the call and the put option last week, we decided to continue to investigate option or derivative trading. This time we looked at how to program binomial trees in MatLab. For a quick refresher a call option is when the owner of the asset, which can be a stock or portfolio, has the right to buy the asset at the pre-determined strike price. A put option is when the owner has the right to sell the asset at the pre-determined strike price. This means that if the asset price goes below the strike price, the put option is advantageous and makes the owner money. If the asset price goes above the strike price then the call option becomes the money making option.

So what exactly does a binomial tree look like. Below is the an example of a binomial tree when n=3 which means that past n=0 there are three distinct division. In this case at n=0 the price of the asset is 100 dollars. Binomial trees are interesting because unlike the Black-Scholes method, which was previously discussed in last weeks post, the binomial tree provides the user with a flexible way to price the option. In this example at each node the price of the asset has two alternatives it can either go up in price or down in price. Unlike the Black-Scholes model the binomial tree provides a period to period transparency in the price of the option. It is important to note that this binomial tree does not follow the same definition that a typical discrete mathematics textbook would describe a tree. This is because the binomial option pricing tree contains cycles, where as in a strict mathematically definition trees cannot contain cycles.

 (Each number is the possible price of the asset).

The binomial tree is more advantageous in pricing American options. This is because for an American option the owner at any point can exercise the option before the maturity date (when the option expires). The binomial tree works backwards to find the value of the option at T=0, when time is equal to zero. This helps determine how much one is willing to pay for the option. The binomial tree finds the difference in price between the original asset price, and the price at another discrete time interval, using the function max(K-P,0) , where K is the strike price, and P is the price of the asset at T=0. This function max(K-P,0) is used for a put option, and described more in our code.

It is important to discuss the variables the binomial tree is dependent on. It is dependent on $\sigma$ which is volatility, simply it is just the standard deviation in returns of an asset or portfolio, $r$ is risk neutral interest rate, $u$ and $d$ are percentages the asset either goes up or down at each stage of the tree. The formula for backwards induction to find the original value of the asset is…

(1)  

where $V_u$ is the upper price at a discrete stage, and $V_d$ is the lower price at a discrete stage. A discrete stage meaning n=1,2… We know need to look at what p is, well p simulates the geometric Brownian motion of the underlying stock or asset with parameters $r$ and $\sigma$. p is given by the equation below:

(2)  

We have finished the pricing of an American put option for a binomial tree when n=3. And plan on continuing to work on the code for the option pricing method to work for an arbitrary number for n, and hope to have completed by our next post.

The past week has been spent turning the basic ferromagnetic spin model (isingtest.m) created last week into a model of a basic financial system. In the ferromagnetic model, the spins of constituent atoms were represented by discrete spin values of 1 or -1. In this model, we will use the analogous atoms to represent agents in a market. This market will be based on any arbitrary good or stock, and like the ferromagnetic model, adjacent “agents” will be able to share information and influence each other. Here, the tendency for atoms to align with adjacent atoms represents a tendency of financial agents that work close to each other to use the same ideas.

In this model, a spin value of 1 signifies that an agent is currently buying the good or stock, whereas the value of -1 signifies that the agent is currently selling it. The market is governed by two rules:

1. The first rule is analogous to the ferromagnetic model: adjacent actors influence each other strongly.
2. The second rule introduces an idea analogous to a heat bath. This rule states that the agents in the minority group have a special kind of knowledge, so all agents in the market will want to be part of that minority group.

In order to initialize this market, we create a grid to store whether or not each agent is a buyer or seller. We initially assigned the majority of agents to be buyers, represented as white pixels. The sellers are interspersed as blocks throughout the market, represented by the black pixels. Below is an example of what an initial market might look like:

From here, our monte carlo method cycles through every agent in the market, using its position relative to its neighbors as well as an overall market trend to calculate the “field” that pushes it toward one strategy or another. After every iteration, the display of the grid was altered accordingly. The equation for this field is below:

(3)  

The first term in the expression is analogous to the spin alignments, whereas the second term is representative to the individual agent’s view of the whole market. The value C corresponds to whether an agent will be a “fundamentalist(C=1)” or a “chartist(C=-1).” If the agent is a fundamentalist, it means that he or she will be likely to follow the rule that it is best to try to get into the minority, making the second term negative and pushing it against the first term. If the agent is a “chartist,” it means that he or she is content with being in the majority, and the sign of the second term is the same as the first. Whether or not the agent is a fundamentalist or a chartist is determined by whether or not choosing this strategy will increase the energy, or risk, of the system, which would improve the likelihood of making money. This probability of an agent changing his or her identity is determined with:

(4)  

Then, the probability of an agent remaining/becoming a buyer is given by:

(5)  

(6)  

The second term is the probability that the agent will remain/become a seller. A random number is generated by the monte carlo method and compared to the probability, deciding which strategy the agent in question will use at the current time. Once the code has cycled through every agent in the market, it will update the grid displaying the identity of each agent’s strategy at the given time. The average value of all of the strategies in the market is then taken, as it is representative of the price of the good or the stock. A positive average represents a positive trend in price, because there are more people that want to buy than sell. The opposite is true as well, as a negative average means that there is a surplus of sellers. Taking this average and subtracting from it the average of the last iteration gives use the change in price, or return, created by a shift in the market. An example is displayed below, plotting this average vs. time:

The periods of extreme return in a negative or positive direction correspond to a market configuration where there are large groups of buyers grouped together and large groups of sellers grouped together as shown below:

Conversely, when the returns stay near 0, it is because there is no real trend in the market and it is more turbulent, as such:

Finally, another way to represent this trend is to show the difference between fundamentalists in the market compared with chartists in the market. This difference is displayed below, and it is possible to see that when there is a more extreme difference, there are more extreme changes in the return.

For our final part of this project, we will try to use the autocorrelation function to characterize the volatility of the good over time. We will also try to add another strategy option, as well as possibly exploring a 3 dimensional system. The code for this part is “TwoStrategyMarketModel.m”

Link to Code:

https://drive.google.com/folderview?id=0B01Mp3IqoCvhflF5azV4bHY2V0c3Vks0VXVyQkIxcnR4RmFaNDhVWEQtZGZSR2t6V2VWOHc&usp=sharing

The files are under the name: BinomialTreeNThree.m and TwoMarketStrategyModel.m

Here is the code for our modeled guitar string and the resulting force on the bridge.  This force oscillates with time, similar to the frequency of the sound waves that the instrument produces.  We were able to find the base and overtone frequencies using the fast fourier transform function in MatLab, however our scaling is off for graphing purposes.  However, we know our data is correct because we plucked the string at a distance (B=1/5) of the total length of the string which resulted in our fourier transform missing a 5th peak.  For future testing, we hope to record a physical string with similar variables to the one in our code, import the sound file to MatLab, analyze the frequencies, and compare it to the results of the modeled string.

——————————————————————————————————————————————

Link to code:

https://docs.google.com/document/d/1vNhrcuUVjmwfNM8raGKmyKAl0Z19bFR8c7rikb9tI1w/edit?usp=sharing

close all
clear all
%%%GUITAR STRING%%%

%% Initialize Variables
L=.65; %Length of string (m)
T=149; %Tension of string (N)
c=320; %Velocity of wave (m/s)
dx=.065; %Distance step
dt=0.00001015625; %time step (s)
r=1; %constant in oscillation calculation equation
Amp=.000005; %Amplitude of pluck (m)
Lpluck=.13; %Length of pluck spot (m)
B=Lpluck/L; %pluck position on string (m)
M=[0]; %Array for y vaules
runtime=1000;

%% Builds initial pluck array
for i=dx:dx:Lpluck
y=i*Amp/Lpluck;
M=[M,y];
end

for i=(Lpluck+dx):dx:L
y=(i-L)*(-Amp/(L-Lpluck));
M=[M,y];
end

%% Build X-aixs
x_axis=[];

for i=0:dx:L
x_axis=[x_axis,i];
end

%% Initialize Loop variables

ynew=zeros(1,length(M));
ycurrent=zeros(1,length(M));
yold=zeros(1,length(M));
ycurrent=M;
yold=M;

F=[];
N=[];

%% Animation Loop (String Oscillations)

for n=1:runtime

for i=2:length(M)-1
ynew(i)=2*(1-r^2)*ycurrent(i)-yold(i)+(r^2)*(ycurrent(i+1)+ycurrent(i-1));
end

f=T*(ynew(2)-ynew(1))/dx;
F=[F,f];

N=[N,(n*dt)];

yold=ycurrent;
ycurrent=ynew;

subplot(3,1,1)
plot(x_axis,ynew)
axis([0,0.7,-Amp*1.1,Amp*1.1])
title(‘Waves on String with Initial Pluck’);
xlabel(‘X displacement (m)’);
ylabel(‘Y displacement (m)’);

subplot(3,1,2)
plot(N,F)
axis([0,runtime*dt,0,0.01])
title(‘Force of String on Bridge’)
xlabel(‘Time (s)’)
ylabel(‘Force (arbitrary units)’)

pause(0.000001)
end

%% Tension/Frequency

FFT=fft(F);

FTscale=2^nextpow2(runtime);
f=dt/2*linspace(0,1,FTscale/2+1);

subplot(3,1,3)
plot(f,(2*abs(FFT(1:FTscale/2+1))))
title(‘Internal Frequencies’)
xlabel(‘Frequency (Hz)’)
ylabel(‘F(t)’)

%need to fix scale on FT

As of 4/13, I have yet to finish hacking and reassembling all of the equipment in my laboratory to take my own EMG (electromyography) data. As a refresher, an EMG is a means of measuring nerve activity by analyzing voltage changes in muscles either on the surface of the skin or subcutaneously. However, I believe that I should be able to complete the relevant assembly of the EMG kit prior to the end of the project.

That being said, I have completed the initial version of my data analysis code to be implemented in this project, as well as in my thesis. This version uses two data set inputs, a raw EMG signal, and a corresponding time vector. In addition, the code requires the sampling frequency. Beyond this, the experimenter has to set noise and contraction thresholds. The noise threshold serves the dual purpose of compensating for any noise not already dealt with in the Fourier transform, as well as setting the region below which the robot will relax and return to the rest state. The contraction threshold, determined by the relative magnitude of successive peaks to a standard strained peak determined immediately following the transform process.

The two Fourier transforms performed on the raw data, forward and back, allow the code to eliminate the most significant predictable sources of noise, that if left ignored, would render any analysis worthless. The three largest sources accounted for currently are movement artifacts, electrical artifacts, and electric noise. The first two noise signals are generated as a byproduct of the EMG process, mostly existing below 20 Hz. Most electric devices emit strong interference at 60 Hz. Following the forward transform, these frequencies are negated and zeroed, then back transformed to create a rectified, useable signal for later analysis.

The code is meant to reliably extract the peak amplitudes and cycle durations of up to two contractions by binning the raw EMG data into smaller datasets that are then analyzed. The reasons for this are that it makes coding less difficult as well as opening up the possibility of a quasi-continuous computation. By chopping  the data into bins and analyzing them separately, the code can process data quickly after the signals are created, which in a living subject, corresponds to continuous motion of the prosthesis. In the picture below, the initial analysis code is represented as a flowchart.

Additionally, I have found a raw sample EMG data set with a time vector, which if needed, can serve as the data for this project. This was available on a public domain, the link to which is below the figure of a raw EMG signal below and in the corresponding google drive document attached to this submission.

http://physionet.org/physiobank/database/emgdb/

Further work will include improving my code such that it is more independent of the data inputted, more modular, and generally faster/more efficient. Additionally, I have yet to comment on the code, which will be another component of this project. I expect for the final data to have improved my code and run analysis on either data that I have taken, or if I cannot get my experiment assembled and working to a reasonable extent, the sample EMG that was available online.

Relevant files can be found here:

https://drive.google.com/drive/u/0/#folders/0B0K4JjJfDy1ffkR4eEZaMzlxb0xvTjRjWDBIUm9iVzZuRmVyeG1aM0RhWHMxcVB6RXZDUFk

Using MATLAB, I’ve created several functions to visualize scattering from a rigid sphere. I’ve included three pictures below, as well as the MATLAB code (Note, code is still in rough form. Commenting and overall readability shall be improved). In all Images, the sound is incident from the left. In the next post, I’ll go into more detail with regard to the actual mathematics; during the past week my greatest priority has been to simply get the code running. For now, here’s the equation describing the sound field due to plane wave scattering from the rigid sphere:

where jn and hn are the Spherical Bessel and Hankel functions of the first kind, respectively, Pn are the associated Legendre Polynomials, x is simply the x coordinate, k is the wave number, and i is the imaginary unit. Notice that since we have an infinite sum, in practice we must truncate and approximate this solution.

The first–a two dimensional cross section of the pressure incident on the sphere from a plane wave (frequency of ~54 Hz). The red curve is simply for visualization of the sphere; it has no other significance in this plot. Note how the greatest pressure is on the face of the sphere, and the least is on the sides and around on the directly opposite side.

Next is the same graph, but extended to 3D. Note the similar shape, but a little more pronounced as this is for 100 Hz and the previous was for ~54 Hz.

Finally, the intensity of the sound field; this includes both incident and scattered sound energy, and one can therefore see the constructive and destructive interference patterns. This is for a plane wave of ~54 Hz.

CODE (note I’ve removed large comments from the front of the code for the purposes of space)

Note: all code to run these simulations should be present. If you can’t get the code to run, please let me know.

Here is a link to the code in Google Drive:

https://drive.google.com/a/vassar.edu/folderview?id=0B3UNwnkM9fLrfjNyTDJEc0tQeklHTTBueFEtMDRZMUVJWnpvV2p4MFpLbHlpVTJ4QXJCXzg&usp=sharing

Alex Molina & Kadeem Nibbs
Computational 375 – Project

GOAL: Week 1: The setup. Establishing the arrays, the initial variables and conditions.

This week we worked on modeling the long jump via MATLAB since we are interested in exploring realistic projectile motion and human acceleration as it pertains to the human body. For this week, we looked to create a baseline code modeling a point particle with air resistance(drag force) and gravity acting on it.

We want to explore the effects of proper body position on the traveled distance in the long jump, so modeling after Chapter 2: Projectile Motion, we worked to generate a realistic curving point through Earth’s atmosphere. Our goal to set a baseline is complete. For the upcoming weeks, we will work to add air resistance and friction in effect. To make the point particle even more interesting, we will start to add cylindrical and spherical shapes in order to resemble a more humanistic body and interesting projectile.

Important Equations:
Drag Force of Air:
Drag Force of Air in the x direction with the x component of velocity
Drag Force of Air in the y direction with the x component of velocity
Important Links:
MatLAB code and other essential information can be found at:
https://drive.google.com/a/vassar.edu/folderview?id=0B6loGd4o7iESfjNoRnVmSjJDNHlabG9qNEJIRmY1Z1JEeS1QNE9rTjlIY2Vqc1NMTWlwdEk&usp=sharing

References:
Giordano, Nicholas J., and Hisao Nakanishi. “Chapter 2: Realistic Projectile Motion.” Computational Physics. 2nd ed. Upper Saddle River, NJ: Pearson/Prentice Hall, 2006. Print.

Our Resulting Projectile Motion via MatLAB (GIF):
We saved the 40 different jpg files from MatLAB and used a gif generator to make this.  It did however erase the line and motion at each point.

Our Resulting Projectile Motion via MatLAB (Actual):

UPCOMING: Week 2: Run simple trial with one-dimensional object in two-dimensional space with air resistance.  Then second trial with a mass-less two dimensional object in 3-dimensional space.

As we began our project, we realized that the Ising model presented in Chapter 8 of Giordano and Nakanishi quickly introduces many aspects that are not relevant to our neural network model, namely the idea of a heat bath at temperature T.  Our neural network is essentially an Ising model assumed to have a temperature of 0, so much of Chapter 8 is unnecessary for our model.  So we skipped ahead and went straight to simple neural networks, as described in Section 12.3 of the text.

The essential elements of our model are an NxN stored pattern, an interaction matrix J_i,j of size N²x N², and an input pattern of size NxN.  After a pattern is stored (by creating the network’s interaction matrix), an input pattern (usually one similar, but not identical, to the stored pattern) is presented.  This pattern is then changed according to the Monte Carlo method, and, if things worked correctly, the output is the same as the stored pattern.

Stored Pattern:

We created a stored pattern on MATLAB by creating a matrix +1’s and -1’s.  Our first pattern is the letter ‘A’, following the examples in Section 12.3 of Giordano and Nakanishi, on a 10×10 matrix.  As mentioned in our earlier posts, each element in the matrix represents a neuron and if the neuron is active it has a value of 1 and if it is inactive, its value is -1. Our ‘A’ pattern is made up of +/-1s where the +1’s represent the pattern and the -1’s represent the background. When these numbers are displayed in the command window, the A pattern is not easy to see, so we use the imagesc command (along with an inverted gray colormap) to better visualize our patterns.

Element Indexing

MATLAB indexes the elements in a matrix with two numbers: the first index refers to the row of the matrix, and the second index refers to the column.  For our neural networks, we want to be able to access every element in the array (twice, separately) in order to have every neuron interact with every other neuron.  The easiest way to do this is by using a single number index, i, which runs from 1 to  N².  In our code, we use i and j.  To actually access the elements in our matrices, we have to use the normal double MATLAB indices, which we call call y and x.  The index i can be calculated from the normal y and x indices using equation 12.20 from the text.

In our code, we frequently have two sets of nested for loops, with which we are able to loop through the entire pattern twice, once with the index i and once with the index j.  For this reason, we use the indices y_i,x_i and y_j,x_j.

Interaction Energies J(i,j):

The interaction energy matrix describes the interaction between every neuron with every other neuron. So if we start off with a P_stored matrix of size NxN, then the J matrix for it is going to be size N²x N². In order to determine the values for the J matrix we need to use equation 12.16 from the text,

where s_i is the value of neuron i in pattern m, and s_j is the value of neuron j in pattern m.  Thus, the J(i,j) matrix is created by multiplying the values of P_stored(y_i,x_i) with the values of P_stored(y_j,x_j). In order to do this on MATLAB we created two for loops where we held the y(rows) values constant and varied our x (columns) values. By doing this we compare P_stored(1,1) to P_stored(1,2), P_stored(1,3), P_stored(1,4)… Then after finishing with the first row the program goes on to compare P_stored(1,1) with P_stored(2,1), P_stored(2,2)… This continues until P_stored(1,1) is compared with all the other points in the matrix and then it does the same thing for P_stored(1,2), P_stored(1,3), and so on.

Example 1:

Below you will see a P_stored matrix and it’s corresponding J matrix. By using equation 12.20 and the double for loop, you are able to index the values for i and j which allows you to create J(i,j).

  Figure 1                        Figure 2

The interaction energy matrix values will later be used to help determine the amount of energy it takes to activate or inactivate a neuron.

Neural Network Model

After finding our interaction energies, we now have an operational memory for our simple neural network. In order for our neural network to be operational we need to create a pattern within P_stored. In order to do this we created an image of the letter A by manually entering in the 1s to create the letter A. We then created a new matrix called P_in to represent the matrix that we will sweep using the Monte Carlo method. P_in is created by randomly flipping the value of some of the neurons in order to distort our original A pattern. In order to make it easier to visualize, we used the function imagesc to turn our matrix values into an image.

After conducting a Monte Carlo sweep over our distorted A-pattern P_in we should get in return the same pattern as our P_stored matrix. Thus, our neural network remembers being fed the desired pattern we stored earlier. In the figures below, the black squares are represented by +1’s in our matrix and the white squares represent the -1’s.

Figure 3

 Figure 4 Figure 5

Monte Carlo Method

The method by which the input pattern is changed into the output pattern is the Monte Carlo method.  The essence of this method is to calculate, for each neuron in the input pattern, the effective energy of the entire neural network in its current state using Equation 12.14.

If the total energy of the network with respect to the neuron i is greater than 0, then neuron i’s value is flipped (+1 -> -1, or vice versa).  If the total energy of the neural network with respect to neuron i is less than 0, its value is unchanged.

This process is carried out for every neuron i, so that every neuron is given a chance to flip.  Once every neuron has been checked, one sweep of the Monte Carlo method is complete.

For our current code, we are only using one Monte Carlo sweep, but as we progress forward we will use multiple Monte Carlo sweeps to check if the output pattern is stable, and to see if multiple steps are needed in order to fully recall the stored pattern.

Our full script that includes all the parts discussed above, is attached below.

Brian and Tewa Initial Code for Neural Networks Project

Future Plans

Next, we will begin to characterize the properties of our model.  We will be investigating questions how different P_in is from the stored pattern, the maximum number of patterns that can be stored in a network, and how long pattern recall can take.  The work for this post was done exclusively together as a team, and moving forward we will begin to split up a bit more as we tackle independent questions.  As we do this, we will be splitting up our code into smaller chunks, some of which will be re-defined as functions, for easier interaction.

The first simulation used in this project was for the precession of the perihelion of Mercury. This involves calculation of Mercury’s orbit using the force law as predicted by general relativity:

We can use the above equation to calculate position and velocity; I used a second-order Runge-Kutta Method. Equations shown are the same for x and y:

Once we have the proper position of Mercury we can convert said position to an angle theta. Then, plotting time against theta we can do a linear fit, the slope of which will be our desired rate of precession. A simulated orbit for Mercury, using alpha = .008, appears as shown:

which already shows the precession of the perihelion over time. In the upcoming weeks I need to more closely investigate the rate of this procession.

The second part of the project involves the investigation of Kirkwood Gaps, or the odd behaviour of solar satellites at certain orbital radii which “happen” to coincide with resonances of Jupiter’s orbit. This meant creating a multi-body simulation (I chose 3-body), so as to calculate the position and velocities of our planets (here, Earth and Jupiter), which are then used to calculate the position and velocities of our asteroids. Luckily, the mass of the asteroid is negligible as compared to that of our largest and most impactful orbiting body, Jupiter, and so the effect of the asteroids on Jupiter’s orbit can be ignored. The initial radius and velocity values used for the asteroids are found in Table 4.4 in Giordano and Nakanishi. Again the Runge-Kutta method was used to compute values. The equations used to calculate these values for planets are similar to the Runge-Kutta equations above and are shown by others doing similar celestial mechanics problems. Using this method I plotted first the orbits of Earth and Jupiter, to make sure my simulation worked:

and then plotted the asteroids’ orbits:

The blue orbit above is for the asteroid found around the 2/1 resonance orbit (ratio of axes is 2:1 for Jupiter:asteroid); here is a visual representation of the effect Jupiter can have on other objects, significantly disturbing their orbit. The next step in this program is to add complexity- more asteroids and a better graphical representation of where these orbital radii “gaps” occur.

Matlab Code:

https://docs.google.com/document/d/1hHy4D-B7eWIbknoHlLHJKvlyEtecRhIF7wupYeuFF1s/

https://docs.google.com/document/d/1jOSlib_gAtcADT_Mc0qFi2PA8n1I2nPnxqqOk_0f34w/

This week I worked on the modeling of Mercury’s orbital motion around the Sun. Using the Euler-Chromer Method and the general relativistic force law (see previous post for a more quantitative description) I was able to create a code that shows Mercury’s movement over several orbits. The value for alpha was used at 0.01 (very much greater than the true value) and that the initial conditions used were when , and . It is clear over several orbits that the location of Mercury’s perihelion is changing.

In the following plot, the different positions for the aphelion (chosen for clarity) are clearly shown by lines connecting the origin (or the Sun) to the aphelion.

It is my next goal to determine the amount of degrees by which the orientation of this orbit changes over time. I have run into the issue that Matlab’s acos and asin (arccos and arcsin) functions give imaginary values if their inputs are not with the range of -1<x<1. When this issue is resolved, I will be able to finish deriving the precession rate of Mercury, and will proceed to experimenting with the eccentricity of Mercurian planet orbits and their resulting rates of precession.

Here is a link to my Matlab code: https://docs.google.com/a/vassar.edu/document/d/1fw8L5ZvhqMMVMYltkJH6Ac425KrmnbbDUJs_6GV7aPI/edit?usp=sharing

(Note that equations are rendered in LaTeX which I typed; they simply come through as images, though I include code for documentation.)

My preliminary results were fairly boring, in that I simply executed in MATLAB the guitar simulation central to the first part of chapter 11 of the text. However, this problem was not trivial, and the machinery constructed here will prove critical later in the project.

As background, recall in chapter 6 how the motion of a string was simulated using eq. 6.6:

$latex y(i,n+1)=2[1-r^2]y(i,n)-y(i,n-1)+r^2 [y(i+1,n)+y(i-1,n)]$

where c = speed of wave and r = cΔt/Δx. The force on the bridge (defined here to be the left boundary of the string) can be calculated using eq. 11.4:

$latex F_{\text{bridge}}=T\frac{y(2,n)-y(1,n)}{\Delta x}$

where T is tension, and y(1,n) is the leftmost element on the string at time n (this is a slight modification to the equation in the book, as MATLAB starts indexing at 1, whereas the book’s authors start at 0). For a string plucked at one fifth its length away from the bridge, the following results were obtained after simulating the vibrating string, calculating the force on the bridge with time, and taking the power spectrum of the force using the discrete Fourier transform written for class:

Figure 1: Preliminary data. Bridge power spectrum calculated in MATLAB matches the results in the book very well, indicating that the code used is a good starting point for the rest of the project.

Code used attached here via DropBox (no need to sign up for an account).