Welcome to TiddlyWiki created by Jeremy Ruston; Copyright © 2004-2007 Jeremy Ruston, Copyright © 2007-2011 UnaMesa Association
!! Monte Carlo method
* [[The FERMIAC|https://en.wikipedia.org/wiki/FERMIAC]]. Metropolis came up with the name "Monte Carlo" with Stan Ulam. Shown is the image of Fermi's mechanical calculator to compute neutron diffusion using Monte Carlo.
!! Literature
# Wm. G. Hoover and C.G. Hoover. Ergodicity of a ~Time-Reversibly Thermostated Harmonic Oscillator and the 2014 Ian Snook Prize. [[(2014)|http://arxiv.org/abs/1408.0256]]. Shuichi Nose opened up a new world of atomistic simulation in 1984. He formulated a Hamiltonian tailored to generate Gibbs' canonical distribution dynamically. This clever idea bridged the gap between microcanonical molecular dynamics and canonical statistical mechanics.
# Peter Olver and Ari Stern. Dispersive Fractalization in linear and nonlinear Fermi–Pasta–Ulam–Tsingou lattices [[arXIv (2020)|https://arxiv.org/abs/2005.12260]].
* The status quo for learning C++: http://www.learncpp.com/
* Great option for those who are confident with other programming languages but need to transfer that to knowledge C++: http://www.cplusplus.com/doc/tutorial/
* https://www.sololearn.com/Course/CPlusPlus/ uses a very interactive learning structure with short readings and examples followed by questions which one must get right to move to the next part of the tutorial.
* http://www.cprogramming.com/tutorial/c++-tutorial.html.
* http://www.cprogramming.com/tutorial/c-tutorial.html
* http://www.learn-c.org/
The innards of a computer begin in its motherboard. This is the base where the CPU, memory, chipset, network, and connectivity slots are located. An annotated image can be found [[here|https://www.getdroidtips.com/types-motherboard-complete-guide/]].
Photo of Intel's "Next Unit of Computing" [[(NUC)|https://en.wikipedia.org/wiki/Next_Unit_of_Computing]] motherboard:
[img(80%+,auto)[http://homepages.spa.umn.edu/~vinals/tspot_files/phys4041/2020/int-nuc-de3815tykhe-3.jpg]]
As shown in the motherboard images, the key part of a computer is the so called processor. They are built in large microfabriaction facilities, contain approximately 1 billion transistors, with individual features down to about 14 nm (and shrinking, 5 nm is the latest technology).
A somewhat modern architecture is the Intel Sandy Bridge processor (with 32 nm features). The processor containes 6 cores or processing elements that share the same 'die", its own cache and memory controleler. An image can be found at http://www.hardwarezone.com.sg/product-intel-core-i7-3960x-extreme-edition . Also see [[Intel Core i9 and AMD Ryzen|https://www.adslzone.net/app/uploads/2017/08/amd-threadripper-intel-core-i9-x.jpg]]
Communications with memory, bus to video cards etc. and with other peripherals is handled by the "chipset". This figure shows a typical schematic layout.
| [img(50%+,auto)[https://upload.wikimedia.org/wikipedia/commons/5/51/Chipset_schematic.svg]] | [img(80%+,auto)[http://homepages.spa.umn.edu/~vinals/tspot_files/phys4041/2020/storage_hierarchy.png]] |
|Communications with memory, bus to video cards etc. and with other peripherals is handled by the "chipset". This figure shows a typical schematic layout ([[Image Credit|https://en.wikipedia.org/wiki/Northbridge_(computing)]])|Storage is hierarchical relative to the "distance" (access time from the processor). Fastest access is more expensive, and hence the hierarchy: from 3-4 clock cycles for L1 cache, to about 20 for L2, and about 90 for L3. ([[Image Credit|https://cs61.seas.harvard.edu/site/2018/Storage2/]])|
Data from cache is not read one number at a time, but rather in blocks ("lines"). Lines are transferred to main memory in "pages". This is done to optimize transfer rates and keep the CPU registers busy. Modern CPUs optimize this processes, going as far as anticipating the data that the processor will need next. Peak transfer to ~DDR4 memory at the moment is around 20 GB/s, with a latency of about 10 ns. As comparison, transfer rates to hard drives (SATA 3.0) are 600 MB/s, or 2000 MB/s with the newer SATA 3.0. For most scientific applications, the hardware component that limits the speed of a computation is memory access, not the speed of the processor per se.
More details can be found in: [[A beginner's guide to High Performance Computing|http://www.shodor.org/media/content/petascale/materials/UPModules/beginnersGuideHPC/moduleDocument_pdf.pdf]] by Rubin Landau, Oregon State U.
Pick one problem below to work on during the semester.
Please note that in a few cases, it is possible to find the relevant code somewhere in the web. Use (or copy, or adaptation) of that code is not acceptable. You are expected to develop your own algorithm in conjunction with what has been discussed in class, write your own code, and perform the analysis of the data.
The grade given for this portion of the class will correspond to the level of sophistication and understanding of your particular method of solution, the algorithm and the code that you have developed, and your understanding of the physical results obtained. In addition, you need to be able and prepared to potentially explain your work orally.
The mid term exam will cover a description of the physical problem, defining mathematical equations (if appropriate), algorithm, and any accuracy or stability considerations. ''A written report covering these aspects of the problem is due on October 31, 2020''.
''A final report containing the description of the algorithm, the file(s) containing the code, and numerical results as appropriate are due on December 17, 2020''.
!!!! Solution of Ordinary Differential Equations
# [[A Renormalization Group Approach to Spontaneous Stochasticity|https://arxiv.org/abs/2007.01333]]. Gregory Eyink, 2020. Bead rolling down a bent wire, fluctuations, and the emergence of spontaneous stochasticity. Solve the equations for the bent wire in Section IV.A as in Sec. IV.B. See also Appendix C.
# ~Fermi-Pasta-Ulam Chain. Basic information in [[Wikipedia|https://en.wikipedia.org/wiki/Fermi%E2%80%93Pasta%E2%80%93Ulam%E2%80%93Tsingou_problem]]. Integrate the system of ordinary differential equations on a discrete lattice. Also solve the continuum PDE limit: the Korteweg–de Vries equation, and obtain numerically a soliton solution.
# Feigenbaum transition to chaos in the logistic map. [Landau], chapter 12. Write the code that finds the bifurcation diagram of the map, and describe the chaotic state.
# Lorentz model. Very simple model of weather forecasting that led to interesting chaotic behavior. The three differential equations that defined the model are given in Eq. (12.50) of [Landau]. [[Wikipedia|https://en.wikipedia.org/wiki/Lorenz_system]] also has a summary introduction. Write your own code to integrate the system of equations, and use it to study the properties of the chaotic state.
!!!! Molecular Simulation
# Molecular Dynamics. Write a Molecular Dynamics code to simulate the equilibrium properties of a fluid. The code should simulate a three dimensional system. The method will be discussed in class, and a two dimensional version will be assigned as a homework. The method is also described in [Landau], Chapter 16.
# [[Diffusion Limited Aggregation|http://www.gut-wirtz.de/dla/]], and the formation of fractal clusters. See also article by [[Thomas Halsey|https://physicstoday.scitation.org/doi/10.1063/1.1333284]].
# [[Coarse grained methods of polymer chain simulation|https://www.sciencedirect.com/science/article/pii/S0021999107005049]].
# [[Brownian dynamics of a set of particles|https://aip.scitation.org/doi/abs/10.1063/1.436761]] including hydrodynamics.
# [[Molecular Dynamics simulation|https://arxiv.org/abs/2009.00216]] of continuum chains. Fourier Transform method.
# Molecular Dynamics for a nematic liquid crystal. Use the [[Gay-Berne|https://www.tandfonline.com/doi/full/10.1080/02678290601140456]] potential as a model of interaction potential of liquid crystal molecules.
!!!! Machine Learning approaches
From the collection on [[Machine Learning for Physical Systems|https://www.sciencedirect.com/journal/journal-of-computational-physics/special-issue/10QSWHT6XHF]].
# [[Parameter reduction in complex models|https://reader.elsevier.com/reader/sd/pii/S0021999119302487?token=FE069D5EA404993DA7D2C4DC42AA2113B841E1F62E6B4AED4384252B944CC5A79BA7A94B459AB960A6F250A76A26E960]].
# [[Data driven identification of physical processes|https://reader.elsevier.com/reader/sd/pii/S002199911930333X?token=45EE6FF2E540779FEBCEC9DAFE1267F9C76B5547861A092949865B55CD6DE52E9BE769BD2C925311450B8C0A8AECE9FB]].
# [[Data driven reconstruction of nonlinear dynamics|https://reader.elsevier.com/reader/sd/pii/S0021999119304474?token=D78CE52D52D5081259C7C1C4A472401014D56938CBA4D6A29ADF444D5CF94554A9B58C723C9EBC22122D0598F2486D52]].
# [[Solving the Schroedinger equation with a neural network|https://reader.elsevier.com/reader/sd/pii/S0021999119306345?token=7483EFAB968DC489BE590FBF380C7937B3762F01A983E3434263B24CDC4ACC023A3CA238D79C97C520E6111DCC02B1CE]].
!!!! Monte Carlo Methods
# Weighted Ensemble Monte Carlo simulation method. [[Daniel Zuckerman and collaborators|https://aip.scitation.org/doi/full/10.1063/1.4821167]]. Implement the method on the enzymatic cycle of Eq. (6) of this manuscript.
# [[The Gillespie algorithm to solve networks of coupled reaction equations|https://en.wikipedia.org/wiki/Gillespie_algorithm]].
# Write a Monte Carlo code, similar to that in [[this page|LabSession 13 2018]], but for the ~Swift-Hohenberg equation. Find the stationary distribution. The required energy is the function $\Phi$ in http://scholarpedia.org/article/Swift-Hohenberg_equation .
# Replica Monte Carlo and Parallel Tempering. [[Section 6|https://arxiv.org/abs/0905.1629]]. A different, more complicated [[example|https://www.sciencedirect.com/science/article/pii/S002199911100619X]].
!!!! Fourier analysis
# [[Ewald Method|http://homepages.spa.umn.edu/~vinals/tspot_files/phys4041/2014/EwaldSum.pdf]] to calculate sums over the Brillouin zone. Compute the electrostatic potential for a given charge distribution.
# Fourier analysis and the Fast Fourier Transform (FFT). [Landau], Chapter 10. Sections 10.2, 10.4. Find a Fourier transform routine on the web, and use it to compute derivatives of a given function.
!!!! Solution of Partial Differential Equations
# See (2), part 2 above under "Solving Ordinary Differential Equations".
# Fluid mechanics with moving boundaries. Use the "Level Set Method" to study the motion of [[bubbles in a fluid|https://www.sciencedirect.com/science/article/pii/S0021999111006218]].
# A numerical algorithm to study phase change. Implement this [[numerical algorithm|https://www.sciencedirect.com/science/article/pii/S037843710900079X]] (given in Eq. (8)) in two spatial dimensions. Study the formation of spatial patterns.
# [[Explicit Runge–Kutta method applied to the incompressible Navier–Stokes equation|https://www.sciencedirect.com/science/article/pii/S0021999111006838]].
# @@color(red):New@@ Liquid crystal defects with regularization.
# @@color(red):New@@ Finite Element Method.
# @@color(red):New@@ Simulated annealing.
# @@color(red):New@@ [[Approximation of slow and fast dynamics in multiscale dynamical systems by the linearized Relaxation Redistribution Method|https://www.sciencedirect.com/science/article/pii/S0021999111006498]].
# @@color(red):New@@ Lattice Boltzmann.
# @@color(red):New@@ [[Sensitivity analysis of limit cycle oscillations|https://www.sciencedirect.com/science/article/pii/S0021999112000071]].
Given the data $\vec{x} = \{1.0,1.3,2.2,2.6,2.8,5.0,7.3,7.4,7.5,7.7,7.9\}$, two Gaussian clusters $k=2$, initial guesses for cluster means $\mu_{1} = 6.63$ and $\mu_{2} = 7.57$, guesses for variances $\sigma_{1} = \sigma_{2} = 1.0$ and weights $c_{1} = c_{2} = 0.5$, the following code performs the iteration:
{{{
PROGRAM em_clustering
! Example 13.4, p. 346. Data mining and analysis. Mohamed Zaki and Wagner
! Meira. Cambridge University Press, 2014.
REAL(8) :: x(11),mu(2),sigma(2),c(2)
REAL(8) :: wg(2,11),w(2,11)
DATA x / 1.0,1.3,2.2,2.6,2.8,5.0,7.3,7.4,7.5,7.7,7.9 /
tpi2=sqrt(8.0*atan(1.0))
mu(1)=6.63
mu(2)=7.57
sigma(1)=1.0
sigma(2)=1.0
c(1)=0.5
c(2)=0.5
! Set up the cluster Gaussians for the current values of the parameters
DO i=1,2
DO j=1,11
wg(i,j)=exp(-(x(j)-mu(i))**2/(2.0*sigma(i)**2))/(tpi2*sigma(i))
END DO
END DO
! Use Bayes theorem to compute the probabilities of the parameters given the data
DO i=1,2
DO j=1,11
w(i,j)=wg(i,j)*c(i)/(wg(1,j)*c(1)+wg(2,j)*c(2))
END DO
END DO
! Recompute means, variances and weights with the newly computed probailities of the parameters
DO i=1,2
sumw=0.0
mu(i)=0.0
DO j=1,11
sumw=sumw+w(i,j)
mu(i)=mu(i)+w(i,j)*x(j)
END DO
mu(i)=mu(i)/sumw
END DO
DO i=1,2
sumw=0.0
sigma(i)=0.0
DO j=1,11
sumw=sumw+w(i,j)
sigma(i)=sigma(i)+w(i,j)*(x(j)-mu(i))**2
END DO
sigma(i)=sigma(i)/sumw
END DO
DO i=1,2
c(i)=0.0
DO j=1,11
c(i)=c(i)+w(i,j)
END DO
c(i)=c(i)/11.0
END DO
DO i=1,2
write(6,9999) i, mu(i),sigma(i),c(i)
9999 FORMAT(i4,4x,3f10.3)
END DO
END
}}}
[img(100%+,auto)[http://homepages.spa.umn.edu/~vinals/tspot_files/em_clustering.jpg]]
[img(100%+,auto)[http://homepages.spa.umn.edu/~vinals/tspot_files/em_clustering_2d.jpg]]
[img(100%+,auto)[http://homepages.spa.umn.edu/~vinals/images/PANO_20130902_180940_rotated2.jpg]]
!Computing Environment
As is standard in modern scientific computing, this class assumes that all work is done in a Unix environment, such as the one provided in the [[School of Physics Unix systems|http://zzz.physics.umn.edu/computing/home]], and in the CSE laboratories (below). Unix provides access to compilers for the high level languages used in the class, to external scientific libraries, and to plotting and visualization software. However, many if not all the tasks associated with this class can be carried out on your personal computer assuming that you install the necessary software (language compilers for the most part). While we can provide some limited help in setting up your personal system, such help is really outside the scope of this class.
It is assumed that you are already familiar with a computing environment, including file editors to create code, at least one high level programming language (Fortran, C, C++, or Python), and the necessary command environments to create code, execute it, and analyze the results. Some limited assistance can be provided here as well, but note that PHYS/AST 4041 is not a programming class.
!!! CSE Computer Lab information
* General [[information|http://cselabs.umn.edu/]], including [[account creation|http://cselabs.umn.edu/accounts]].
* ''Remote access'' ([[offsite|http://help.cselabs.umn.edu/offsite]]) to Labs:
** Browser based access to Unix environment: https://vole.cse.umn.edu
** It is possible to connect directly to desktops in CSE Computer Labs. [[Information|https://cseit.umn.edu/]], including [[ssh|https://cseit.umn.edu/knowledge-help/remote-linux-applications-over-ssh]] and [[VPN services|https://cseit.umn.edu/knowledge-help/establish-vpn-connection]]. The list of available individual desktops is [[here|https://cseit.umn.edu/computer-classrooms/cse-labs-unix-machine-listings]].
* ''Local access'' to Unix environment from the actual laboratories.
** At the Start menu (left hand side bottom), start the ~FastX 2 application.
** Click on + sign (right side of menu), and click on web.
** Name: enter any name for session. Host: https://csel-vole-nn.cselabs.umn.edu:3443 nn are the two digits shown on the label of your screen monitor (lower left). User: ~YourUsername. Save.
** Double click on the connection just created. Enter your username and password at the prompts. Click continue on the next window that appears.
** Click on + sign on top right (create a connection). Eventually a "deskpace image" will appear. Double click on it. Now you have a Unix session running. You can maximize the window if you want it to take the whole screen.
** To terminate, just click on window icon x
** At the bottom of the Unix desktop, you will find an icon for a "terminal" and one for a browser. You can click on them to open.
! Programming languages
In order to follow current standards in scientific computing, we recommend that you use at least one of the following high level languages in homework assignments or lab sessions: C, C++ or FORTRAN. We will also use Python in the context of short demonstrations. If you are more comfortable, you are welcome to use MATLAB, Java, etc. for your own development needs. However, by the end of the semester you should have become comfortable with one of the other four.
CSE Labs provide development tools, compiler, and libraries for all four. If you choose to use your personal computer, you will have to download the respective compilers. The compilers for all four languages are open source and free.
!!! C/C++
* Compiler: g++ filename.c or gcc filename.c -lm (type man g++ or man gcc for information and options).
* [[C tutorials|CTutorials]], [[C++ tutorials|C++Tutorials]].
* Analysis of [[common features|http://www.featflow.de/beta/en/software/featflow2/tutorial/tutorial_lang.html]] between FORTRAN and C++.
!!! FORTRAN
* Compiler: gfortran filename.f90 (type man gfortran for information and options).
* [[FORTRAN tutorials|http://fortranwiki.org/fortran/show/Tutorials]]. FORTRAN reference [[here|http://www.am.qub.ac.uk/users/d.dundas/msci_workshop/fortran/notes/fortran90.pdf]] and [[here|http://www.personal.psu.edu/jhm/f90/lectures/]]. See also [[AMath 483|http://faculty.washington.edu/rjl/uwamath583s11/sphinx/notes/html/]] at University of Washington, including an introduction to Unix.
* [[Introduction|http://www.cs.rpi.edu/~szymansk/OOF90/F90_Objects.html]] to object oriented FORTRAN 90.
!!! Python
* [[Introduction|http://faculty.washington.edu/rjl/uwamath583s11/sphinx/notes/html/]].
* [[How to use FORTRAN subroutines in python| https://github.com/lucasmyers97/F2PYInstructions]] (by Lucas Myers).
!!! Other software resources
* [[Anydesk|https://anydesk.com/en/downloads/linux]] or [[Teamviewer|https://www.teamviewer.com/en/]] should be installed on your personal system in order to interact with the instructor during Laboratory Sessions. The software allows remote access to your desktop, and therefore allows the instructor to provide help with your work.
* Data plotting and visualization: [[gnuplot|http://www.gnuplot.info/]] and [[xmgrace|http://plasma-gate.weizmann.ac.il/Grace/doc/UsersGuide.html]]. Both are available on School computers and the CSE Labs.
* MATLAB can also be used for data manipulation and visualization (although do not use MATLAB for your homework assignments). The University has a site license for the package ([[information|http://zzz.physics.umn.edu/computing/software/matlab]]). To use, just type matlab in any School computer.
* Intel Math Kernel Lirbrary, Reference Manual: https://software.intel.com/en-us/mkl_11.2_ref
! ~GitHub
The version control and repository environment ~GitHub will be used to make code and other class materials available to the class. This is also an open source framework that you can install on your personal computer, as well as run it on the CSE Labs server.
* Access: https://github.umn.edu .
* Log in and email the instructor when you have done so. The instructor will then authorize you to access the repository for this class.
* Initialize. Go to (or create) the subdirectory in the CSE Labs computer that you want to use for this class and type
{{{
git clone https://github.umn.edu/vinals/PHYS4041-2020.git
}}}
You can also clone (make a copy) of all the repository content in your home computer or laptop by issuing the same command. Keep in mind that you will need to have git installed on that computer. Versions of git exist for all modern operating systems, including phones (if you care about such things).
To update the local contents on your directory on the CSE Labs host, go to that directory and type
{{{
git pull
}}}
{{justifyright{
''Due Tuesday, September 16, 2014''
}}}
@@color(blue):Important note:@@ Code needs to be submitted by email to the TA. Other material: Hardcopies with plots, calculations, etc. (but ''not code'') are to be submitted directly to the instructor. It is also possible to submit graphs in image form or scripts to generate them directly to the TA in electronic form. It is best if you group all files into a single archive and email the archive. For example, create a directory homework1, copy all files that need to be submitted to it and create the archive by typing:
* % mkdir homework1
* % cp [filename1] [filename2] ... [filenamen] homework1/
* % tar cvf homework1.tar homework1
and email file homework1.tar
<hr>
# Write a subroutine or a function that takes an array of real numbers as argument, and that returns +1 if the sum of the array elements is positive, 0 if the sum of array elements is zero, and -1 if the sum is negative. Write a main program that uses this subroutine or function on a particular example of your choice.
# Consider the following system of two equations with two unknowns \begin{eqnarray}
3x + ay & = & 10, \\
5x + by & = & 20,
\end{eqnarray}
with $a = 2.100 \pm 5 \times 10^{-4}$ and $b = 3.300 \pm 5 \times 10^{-4}$. Solve for $y$ and find the error in $y$.
# Given a sequence of real numbers $\{x_{1}, x_{2}, ... x_{n} \}$, devise and describe an iterative algorithm to sort the sequence in (a) ascending value, (b) ascending absolute value. Write a computer program to implement your algorithm, and run it on a test sequence of your choosing.
# Write some code to establish the number of digits of precision available in both single and double precision floating point arithmetic in the computer that you are using. To accomplish this task, obtain the following results with successively smaller values of $\epsilon$: $x = 1 + \epsilon$, $y = 1 - \epsilon$ and $z = x - y$. Mathematically, $z = 2\epsilon$, but as $\epsilon$ becomes sufficiently small, the code will return $z = 0$ instead. Find for both cases of single and double precision how small can $\epsilon$ be until you obtain $z = 0$; this is the precision of the computer you are using.
# Write some code that uses single precision arithmetic to evaluate the sum
$$ S_{N} = \sum_{n=1}^{2N} (-1)^{n} \frac{n}{n+1} = - \sum_{n=1}^{N} \frac{2n-1}{2n} + \sum_{n=1}^{N} \frac{2n}{2n+1} = \sum_{n=1}^{N} \frac{1}{2n(2n+1)}.
$$
## Evaluate the sum $S_{N}$ by using the three algebraic expressions given above when N = 100, 1000, 10000. Comment on any differences you may find among the three calculations.
## Using the third algebraic formula determine how many terms would you require in the calculation to make sure that the sum has converged to an accuracy of $10^{-7}$.
[[pdf version|http://homepages.spa.umn.edu/~vinals/tspot_files/phys4041/2014/PHYS4041_Homework1.pdf]]
{{justifyright{
''Due Friday, September 23, 2016''
}}}
@@color(blue):Important note:@@ Code needs to be submitted by email to the TA. Other material: Hardcopies with plots, calculations, etc. (but ''not code'') are to be submitted directly to the instructor or TA. It is also possible to submit graphs directly to the TA in electronic form. It is best if you group all files into a single archive and email the archive. For example, create a directory homework1, copy all files that need to be submitted to it and create the archive by typing:
{{{
mkdir homework1
cp [filename1] [filename2] ... [filenamen] homework1/
tar cvf homework1.tar homework1
}}}
and email file homework1.tar
<hr>
# Write a subroutine or a function that takes an array of real numbers as argument, and that returns +1 if the sum of the array elements is positive, 0 if the sum of array elements is zero, and -1 if the sum is negative. Write a main program that uses this subroutine or function on a particular example of your choice.
# Consider the following system of two equations with two unknowns \begin{eqnarray}
3x + ay & = & 10, \\
5x + by & = & 20,
\end{eqnarray}
with $a = 2.100 \pm 5 \times 10^{-4}$ and $b = 3.300 \pm 5 \times 10^{-4}$. Solve for $y$ and find the error in $y$.
# Given a sequence of real numbers $\{x_{1}, x_{2}, ... x_{n} \}$, devise and describe an iterative algorithm to sort the sequence in (a) ascending value, (b) ascending absolute value. Write a computer program to implement your algorithm, and run it on a test sequence of your choosing.
# Write some code to establish the number of digits of precision available in both single and double precision floating point arithmetic in the computer that you are using. To accomplish this task, obtain the following results with successively smaller values of $\epsilon$: $x = 1 + \epsilon$, $y = 1 - \epsilon$ and $z = x - y$. Mathematically, $z = 2\epsilon$, but as $\epsilon$ becomes sufficiently small, the code will return $z = 0$ instead. Find for both cases of single and double precision how small can $\epsilon$ be until you obtain $z = 0$; this is the precision of the computer you are using.
# Write some code that uses single precision arithmetic to evaluate the sum
$$ S_{N} = \sum_{n=1}^{2N} (-1)^{n} \frac{n}{n+1} = - \sum_{n=1}^{N} \frac{2n-1}{2n} + \sum_{n=1}^{N} \frac{2n}{2n+1} = \sum_{n=1}^{N} \frac{1}{2n(2n+1)}.
$$
## Evaluate the sum $S_{N}$ by using the three algebraic expressions given above when N = 100, 1000, 10000. Comment on any differences you may find among the three calculations.
## Using the third algebraic formula determine how many terms would you require in the calculation to make sure that the sum has converged to an accuracy of $10^{-7}$.
[[pdf version|http://homepages.spa.umn.edu/~vinals/tspot_files/phys4041/2016/PHYS4041_Homework1.pdf]]
{{justifyright{
''Due Tuesday, September 18, 2018''
}}}
@@color(blue):Important note:@@ Code needs to be submitted by email to the TA. Other material: Hardcopies with plots, calculations, etc. (but ''not code'') are to be submitted directly to the instructor or TA. It is also possible to submit graphs directly to the TA in electronic form. It is best if you group all files into a single archive and email the archive. For example, create a directory homework1, copy all files that need to be submitted to it and create the archive by typing:
{{{
mkdir homework1
cp [filename1] [filename2] ... [filenamen] homework1/
tar cvf homework1.tar homework1
}}}
and email file homework1.tar
<hr>
!!! Question 1. Attainable accuracy in single and double precision arithmetic
Write and use code to establish the number of available digits in both single and double precision floating point arithmetic in the computer that you are using (state what computer it is). In order to do this, write a code that computes the following quantities as a function of $\epsilon$: $x= 1.0 + \epsilon$, $y = 1.0 - \epsilon$, $z = x -y$. The exact result is $z = 2 \epsilon$. However, as you run your code for smaller and smaller values of $\epsilon$, it will eventually return $z = 0$. At the point, the value of $\epsilon$ used is the precision of your computer.
!!! Question 2. Practice function coding
Write a subroutine or a function that takes an array of real numbers as argument, and that returns +1 if the sum of the array elements is positive, 0 if the sum of array elements is zero, and -1 if the sum is negative. Write a main program that uses this subroutine or function on a particular example of your choice.
!!! Question 3. Error propagation
Consider the following system of two equations with two unknowns \begin{eqnarray}
3x + ay & = & 10, \\
5x + by & = & 20,
\end{eqnarray}
with $a = 2.100 \pm 5 \times 10^{-4}$ and $b = 3.300 \pm 5 \times 10^{-4}$. Solve for $y$ and find the error in $y$.
!!! Question 4. Error accumulation in series evaluation
Write some code that uses single precision arithmetic to evaluate the sum
$$ S_{N} = \sum_{n=1}^{2N} (-1)^{n} \frac{n}{n+1} = - \sum_{n=1}^{N} \frac{2n-1}{2n} + \sum_{n=1}^{N} \frac{2n}{2n+1} = \sum_{n=1}^{N} \frac{1}{2n(2n+1)}.
$$
## Evaluate the sum $S_{N}$ by using the three algebraic expressions given above when N = 100, 1000, 10000. Comment on any differences you may find among the three calculations.
## Using the third algebraic formula determine how many terms would you require in the calculation to make sure that the sum has converged to an accuracy of $10^{-7}$.
!!! Question 5. Algorithm to conduct an iteration
Given a sequence of real numbers $\{x_{1}, x_{2}, ... x_{n} \}$, devise and describe an iterative algorithm to sort the sequence in (a) ascending value, (b) ascending absolute value. Write a computer program to implement your algorithm, and run it on a test sequence of your choosing.
{{justifyright{
''Due Tuesday, September 22, 2020''
}}}
[[pdf version|http://homepages.spa.umn.edu/~vinals/tspot_files/phys4041/2020/hwk1.pdf]]
@@color(blue):Important note:@@ Use canvas to upload your homework. If it is hand written work, please scan and upload a .pdf file. If it is code, you can simply upload the file(s).
<hr>
!!! Question 1. Practice function coding
Write a subroutine or function that takes an array of real numbers as argument, and then returns +1 if the sum of the array elements is positive, 0 if the sum of array elements is zero, and -1 if the sum is negative. Write a main program that uses this subroutine or function on a particular example of your choice.
!!! Question 2. Error propagation. Pen and paper !
Consider the following system of two equations with two unknowns \begin{eqnarray}
3x + ay & = & 10, \\
5x + by & = & 20,
\end{eqnarray}
with $a = 2.100 \pm 5 \times 10^{-4}$ and $b = 3.300 \pm 5 \times 10^{-4}$. Analytically solve for $y$, and find the error in $y$.
!!! Question 3. Error accumulation in series evaluation
Write some code that uses single precision arithmetic to evaluate the sum
$$ S_{N} = \sum_{n=1}^{2N} (-1)^{n} \frac{n}{n+1} = - \sum_{n=1}^{N} \frac{2n-1}{2n} + \sum_{n=1}^{N} \frac{2n}{2n+1} = \sum_{n=1}^{N} \frac{1}{2n(2n+1)}.
$$
## Evaluate the sum $S_{N}$ by using the three algebraic expressions given above when N = 100, 1000, 10000. Comment on any differences you may find among the three calculations.
## Using the third algebraic formula determine how many terms would you require in the calculation to make sure that the sum has converged to an accuracy of $10^{-7}$.
!!! Question 4. Algorithm to conduct an iteration
Given a sequence of real numbers $\{x_{1}, x_{2}, ... x_{n} \}$, devise and describe an iterative algorithm to sort the sequence in (a) ascending value, (b) ascending absolute value. Write a computer program to implement your algorithm, and run it on a test sequence of your choosing.
!!! Question 5. Select term project
Go to [[Exam Projects]] and select your term project. Consult or discus with the instructor as needed.
{{justifyright{
''Due Friday, December 11, 2018''
}}}
!!Monte Carlo simulation
Complete the development of a Monte Carlo algorithm to simulate the equilibrium distribution of the Ising model in two dimensions just as done in class. Consider a square lattice with $N \times N$ sites, in which each lattice site $(i,j)$ is occupied by a spin variable adopting the values $S{ij} = \pm 1$.
* The energy of a given configuration $\{ S_{ij} \}$ is given by, $$ E(\{ S_{ij} \}) = -J \sum_{<n.n.>} S_{ij}S_{kl}$$
where the constant $J$ is the interaction constant, and the sum extends only over nearest neighbors in the square lattice, without double counting interection bonds (there are four nearest neighbors in the square lattice: north, south, east and west of the spin in question).
* Use dimensionless variables: $J^{\prime} = J/k_{B}T$.
* Use periodic boundary conditions in both spatial directions. The following modular arithmetic may be useful. If $i$ is a given index (in either direction), then the next index, $ip1$ (i plus one), and the prior index $im1$ (i minus 1) can be expressed as:$$ip1 = 1 + MOD(i,N) \quad\quad im1 = 1+MOD(i+N-2,N),$$ where $MOD$ is the MOD function (MOD(a,b) in FORTRAN, a%b in C).
* Starting from a random initial condition, run two simulations for $J^{\prime} = 0.4, 0.5$. Use the Metropolis algorithm as discussed in class.
* Compute the average energy per spin $\langle E \rangle/N^{2}$ and the total magnetization per spin $\langle M \rangle = (1/N^{2}) \langle \sum_{i,j} S_{i,j} \rangle$.
!! Genetic algorithm
Write the code that uses a genetic algorithm to find the maximum of the function $f(x) = x^{2}$ for $x$ an integer in the interval $[0,31]$. Use judicious choices for the sizes of the sample, and the probabilities or frequencies of crossovers and mutations.
{{justifyright{
''Due Tuesday, September 23, 2014''
}}}
!!!!1. Bessel function approximation
Bessel functions are trancendental functions that can be obtained by summing a series. However, there are recursion formulae that make the computation more efficient, although, as seen in class, recursion can introduce its own numerical problems. Use the following recursion formulae to compute the Bessel function $J_{l}$ with $l = 0, 1, 2, 3, 4$ and plot the results for $x \in (0,12)$. The recursion formulae ("upwards" or "downwards") are \begin{eqnarray}
J_{l+1}(x) & = & \frac{2l+1}{x} J_{l}(x) - J_{l-1}(x) \\
J_{l-1}(x) & = & \frac{2l+1}{x} J_{l}(x) - J_{l+1}(x) \\
J_{0}(x) = \frac{\sin x}{x} & \; & J_{1}(x) = \frac{\sin x - x \cos x}{x^{2}}
\end{eqnarray} Using the following table to help evaluate the accuracy of your approximations,
|! x |! $J_{3}$ |! $J_{5}$ |! $J_{8}$ |
| 0.10 | 9. 518519719 10$^{-6}$ | 9.616310231 10$^{-10}$ | 2.901200102 10$^{-16}$ |
| 1.00 | 9.006581118 10$^{-3}$ | 9.256115862 10$^{-5}$ | 2.826498802 10$^{-8}$ |
| 10.0 | -3.949584498 10$^{-2}$ | -5.553451162 10$^{-1}$ | 1.255780236 10$^{-1}$ |
** Compute the Bessel functions required by using the upward recursion.
** Compare your solution to the table. Argue how errors might have propagated and what effect this may have on the accuracy of the solution.
** Chose an arbitrary starting function $J_{L}(x)$ with large $L$ and use the iteration downward. Compare the results to the upward iteration for the Bessel functions asked above.
!!!!2. Spline fitting
We have found experimentally the pressure $p$ as a function of the molar volume ($v$) of a binary mixture at constant temperature and composition (in some arbitrary units). The results are shown in the following table.
|v | 1.0 | 1.1 | 1.2 | 1.3 | 1.4 | 1.5 |
| p | 4.7 | 3.5 | 3.0 | 2.7 | 3.0 | 2.4 |
## Write your own program to find the cubic spline interpolation to $p=p(v)$ in this interval. The description of the algorithm can be found at http://mathworld.wolfram.com/CubicSpline.html
## Use the spline interpolation to find the value of the volume at which $p=3.25$.
## Accomplish the same task but download publicly available code on the web (scroll to the bottom of the page):
** Fortran version: http://people.sc.fsu.edu/~jburkardt/f_src/spline/spline.html
** C version: http://people.sc.fsu.edu/~jburkardt/c_src/spline/spline.html
** C++ version: http://people.sc.fsu.edu/~jburkardt/cpp_src/spline/spline.html
[[pdf version|http://homepages.spa.umn.edu/~vinals/tspot_files/phys4041/2014/PHYS4041_Homework2.pdf]]
{{justifyright{
''Due Friday, September 30, 2016''
}}}
!!! 1. Solution of simple linear system of equations (LAPACK)
Write a code (and makefile) that using the Lapack routine dgesv solves the following system of equations:
\begin{eqnarray}
3 x_{1} + 2 x_{2} & = & 2 \\
2 x_{1} + 6 x_{2} & = & -8
\end{eqnarray}
Documentation for the library can be found [[here|http://www.netlib.org/lapack/double/dgesv.f]]. This document assumes FORTRAN language. If you use C, call dgesv_(&n,&nrhs,a,&lda,ipiv,b,&ldb,&info) (note the underscore). The arguments and usage information are the same in both languages.
!!! 2. Scaling analysis of LAPACK routine dgesv
We wish to examine how the execution time grows with $n$, the order of the system of equations to solve. In order to do so,
# Create a large matrix $A$ with random entries. CALL srand(int seed) (seed is an arbitrary integer) will initialize the random sequence, and rand() returns real random numbers uniformly distributed between 0 and 1. Create a matrix of size $(n,n)$ with random entries, and also a right hand side vector $b$. To make the problem better conditioned, add to the diagonal elemens of $A$ some positive quantity (e.g. 1).
# Solve the linear system of equations for increasing values of $n$. Use the time command to obtain the execution time. For example, if a.out is the executable, type instead
{{{
time a.out
}}}
This will give you the execution time. Plot in a log log scale the execution time with $n$ and examine how the time scales with $n$ for large $n$.
Although simple, this method will also include the time involved in the other tasks in your code. For example, it will include the time spent setting up the matrix A. A more accurate timing can be obtained by using the function cpu_time (in fortran):
{{{
CALL cpu_time(t0)
! ... tasks to be timed ...
CALL cpu_time(t1)
WRITE(6,*) 'Elapsed time (in seconds): ',t1-t0
}}}
In C you would [[use|https://www.gnu.org/software/libc/manual/html_node/CPU-Time.html]]:
{{{
#include <time.h>
clock_t start, end;
double cpu_time_used;
start = clock();
… /* Do the work. */
end = clock();
cpu_time_used = ((double) (end - start)) / CLOCKS_PER_SEC;
}}}
!!! 3. Zero of function with spline interpolation
We have found experimentally the pressure $p$ as a function of the molar volume $v$ of a binary mixture at constant temperature and composition (in some arbitrary units). The results are shown in the following table.
|v | 1.0 | 1.1 | 1.2 | 1.3 | 1.4 | 1.5 |
| p | 4.7 | 3.5 | 3.0 | 2.7 | 3.0 | 2.4 |
# Write your own code to find a cubic spline interpolation polynomial $p=p(v)$ of this data in the interval given. Plot the data along with the spline interpolation on the same graph.
# Use the spline interpolation to find the value of the volume at which $p=3.25$. Do it approximately on the graph first.
# Compute the value of the volume at which $p=3.25$ by solving the equation $f(v)=p(v)-3.25=0$ by using the ~Newton-Raphson method ''WITH'' the interpolated spline function.
# Additional information
## It is not necessary for the solution of this problem, but there is [[documentation|http://mathworld.wolfram.com/CubicSpline.html]] about the definition and calculation of spline interpolating polynomials.
## If you use FORTRAN, you can use the routines spline.f and splint.f provided in the class github repository.
## If you use C, you can use the routines spline.c and splint.c and header files nrutil.h, nrutils.c provided in the class github repository.
{{justifyright{
''Due Tuesday, September 25, 2018''
}}}
!!! 1. Bessel function approximation
Bessel functions are transcendental functions that can be computed by summing a series. However, there are recursion formulae that make the computation more efficient, although, as seen in class, recursion can introduce its own numerical problems. Use the following recursion formulae to compute the Bessel function $J_{l}$ with $l = 0, 1, \ldots, 8$ and plot the results for $x \in (0,12)$. The recursion formulae ("upwards" or "downwards") are \begin{eqnarray}
J_{l+1}(x) & = & \frac{2l+1}{x} J_{l}(x) - J_{l-1}(x) \\
J_{l-1}(x) & = & \frac{2l+1}{x} J_{l}(x) - J_{l+1}(x) \\
{\rm Note:} \;\; J_{0}(x) = \frac{\sin x}{x} & \; & J_{1}(x) = \frac{\sin x - x \cos x}{x^{2}}
\end{eqnarray} The following table lists accurate values of selected Bessel functions, and can be used to help evaluate the accuracy of your approximations,
|! x |! $J_{3}$ |! $J_{5}$ |! $J_{8}$ |
| 0.10 | 9. 518519719 10$^{-6}$ | 9.616310231 10$^{-10}$ | 2.901200102 10$^{-16}$ |
| 1.00 | 9.006581118 10$^{-3}$ | 9.256115862 10$^{-5}$ | 2.826498802 10$^{-8}$ |
| 10.0 | -3.949584498 10$^{-2}$ | -5.553451162 10$^{-1}$ | 1.255780236 10$^{-1}$ |
# Write a code that computes the Bessel functions required by using the upward recursion.
# Compare your solution to the table. Argue how errors might have propagated and what effect this may have on the accuracy of the solution.
# Chose arbitrary starting values for $J_{L}(x)$ and $J_{L-1}(x)$ with large $L$ and use the iteration downward. Then use the ratio between the analytical and recursive values for $J_{0}(x)$ to normalize the value for the desired $J_{l } (x)$. Compare the results to the upward iteration for the Bessel functions asked above.
!!! 2. Solution of simple linear system of equations (LAPACK)
Write a code (and makefile) that using the Lapack routine dgesv solves the following system of equations:
\begin{eqnarray}
3 x_{1} + 2 x_{2} & = & 2 \\
2 x_{1} + 6 x_{2} & = & -8
\end{eqnarray}
Documentation for the library can be found [[here|http://www.netlib.org/lapack/double/dgesv.f]]. This document assumes FORTRAN language. If you use C/C++, call dgesv_(&n,&nrhs,a,&lda,ipiv,b,&ldb,&info) where the length of ''array'' a is $lda*n$ (note also the underscore). Documentation [[here|http://www.physics.utah.edu/~detar/phys6720/handouts/lapack.html]]. The arguments and usage information are the same in both languages. In C++ keep in mind the extra statement needed to indicate that dgesv has been compiled in C.
!!! 3. Scaling analysis of LAPACK routine dgesv
We wish to examine how the execution time grows with $n$, the order of the system of equations to solve. In order to do so, write a program that
# Creates a large real matrix of order $n$, $A(n,n)$, with random entries. In FORTRAN, use CALL srand(seed) (where seed is an arbitrary integer) to initialize the random sequence. In C, use the function void srand(unsigned int seed) instead. Then, in both languages, the function rand() returns real random numbers uniformly distributed between 0 and 1. Create a matrix of size $(n,n)$ with random entries, and also a right hand side vector $b$. To make the problem better conditioned, add to the diagonal elements of $A$ some positive quantity (e.g. 1).
# Solve the linear system of equations for increasing values of $n$. Use the time command to obtain the execution time. For example, if a.out is the executable, type instead
{{{
time a.out
}}}
This will give you the execution time. Plot in a log log scale the execution time with $n$ and examine how the time scales with $n$ for large $n$.
Although simple, this method will also include the time involved in the other tasks in your code. For example, it will include the time spent setting up the matrix A. A more accurate timing can be obtained by using the function cpu_time (in fortran):
{{{
CALL cpu_time(t0)
! ... tasks to be timed ...
CALL cpu_time(t1)
WRITE(6,*) 'Elapsed time (in seconds): ',t1-t0
}}}
In C you would [[use|https://www.gnu.org/software/libc/manual/html_node/CPU-Time.html]]:
{{{
#include <time.h>
clock_t start, end;
double cpu_time_used;
start = clock();
… /* Do the work. */
end = clock();
cpu_time_used = ((double) (end - start)) / CLOCKS_PER_SEC;
}}}
{{justifyright{
''Due Wednesday, September 30, 2020''
}}}
!!! 0. Format requirements
In order to streamline uploading into Canvas and grading, please adhere to the following:
# Submit one separate source file for each of the three questions. If there is a main() and functions for a question, please include them in the same file.
# If you submit a Jupyter notebook, please also submit a .pdf export to simplify quick referencing.
# Results, plots, and discussion should be contained in a separate, single file (any format that you wish. If the format is obscure, please export to .pdf and submit the pdf file).
# Do not submit executables (after compilation). They are not portable in principle, and cannot be used. Include makefiles or specify compilation process where necessary (e.g. when multiple source files are involved).
!!! 1. Scaling analysis of LAPACK routine dgesv
During Lab Session 1, you learned to link and use a LAPACK routine to solve a linear system of equations. We wish now to examine how the execution time used by the routine grows with $n$, the number of equations in the system to solve. In order to do so,
'' Create a large matrix $A$ with random entries''
* In Fortran, {{{CALL srand(seed)}}} (seed is an arbitrary integer) will initialize the random sequence, and {{{rand()}}} returns real random numbers uniformly distributed between 0 and 1.
* In C/C++, {{{rand()}}} returns an integer in the range from 0 to {{{RAND_MAX}}}. Adjust accordingly.
{{{
#include <stdlib.h>
unsigned seed;
srand(seed);
r=rand();
}}}
* In Python, for random floating-point numbers between 0 and 1,
{{{
import numpy as np
seed_value = ...
np.random.seed(seed_value)
r = np.random.rand()
}}}
Create a matrix of size $(n,n)$ with random entries, and also a right hand side vector $b$. To make the problem better conditioned, add to the diagonal elements of $A$ some positive quantity (e.g. 1).
''Solve the linear system of equations for increasing values of $n$ and time the execution''
//In Unix// the {{{time}}} command can be used to obtain the execution time. For example, if a.out is the executable, enter
{{{
time a.out
}}}
This will give you the execution time. Plot in a log log scale the execution time with $n$ and examine how the time scales with $n$ for large $n$.
Although simple, this method will also include the time involved in the other tasks in your code. For example, it will include the time spent setting up the matrix A. A more accurate timing can be obtained by using the function cpu_time (in fortran):
{{{
CALL cpu_time(t0)
! ... tasks to be timed ...
CALL cpu_time(t1)
WRITE(6,*) 'Elapsed time (in seconds): ',t1-t0
}}}
In C you would [[use|https://www.gnu.org/software/libc/manual/html_node/CPU-Time.html]]:
{{{
#include <time.h>
clock_t start, end;
double cpu_time_used;
start = clock();
/* ... Do the work ... */
end = clock();
cpu_time_used = ((double) (end - start)) / CLOCKS_PER_SEC;
}}}
In Python, a portable approach would be
{{{
from timeit import default_timer
start = default_timer();
# ... tasks to be timed ...
end = default_timer();
seconds_elapsed = end - start
}}}
You can also use the [[timeit module|https://docs.python.org/3/library/timeit.html#examples]] to time multiple iterations of the task in question in one go for better accuracy.
!!! 2. Spline Interpolation
We have found experimentally the pressure $p$ as a function of the molar volume ($v$) of a binary mixture at constant temperature and composition (in some arbitrary units). The results are shown in the following table.
|! v | 1.0 | 1.1 | 1.2 | 1.3 | 1.4 | 1.5 |
|! p | 4.7 | 3.5 | 3.0 | 2.7 | 3.0 | 2.4 |
## Write your own program to find the cubic spline interpolation to $p=p(v)$ in this interval. The description of the algorithm can be found at http://mathworld.wolfram.com/CubicSpline.html
## Use the spline interpolation code to find the value of the volume at which $p=3.25$. You will need a routine to solve a tri-diagonal system of equations. You will find it in the class repository under tridag.f, tridag.c, and tridag.py.
!!! 3. Spline interpolation using external libraries
Accomplish the same task but by downloading publicly available libraries on the web or on the class repository (remember to issue git pull, and look under homework 2).
** Fortran version: spline.f, splint.f or http://people.sc.fsu.edu/~jburkardt/f_src/spline/spline.html (scroll to the bottom of the page)
** C version: splint.c, spline.c, and nrutil.h, or http://people.sc.fsu.edu/~jburkardt/c_src/spline/spline.html (scroll to the bottom of the page)
** C++ version: You can use the C routines above, or http://people.sc.fsu.edu/~jburkardt/cpp_src/spline/spline.html (scroll to the bottom of the page)
** Python function: https://docs.scipy.org/doc/scipy/reference/tutorial/interpolate.html#spline-interpolation
{{justifyright{
''Due Tuesday, October 8, 2014''
}}}
!!!!1. Solution of Poisson equation
Poisson's equation relates the electrostatic potential $V(x)$ to a charge distribution $\rho(x)$. In arbitrary units, and in two spatial dimensions we write,
\begin{equation}
\frac{\partial^{2} V(x,y)}{\partial x^{2}} + \frac{\partial^{2} V(x,y)}{\partial y^{2}} = - \rho(x,y)
\end{equation}
By a specific discretization scheme, the equation is written as ($\Delta x = \Delta y$):
\begin{equation}
\frac{1}{\Delta x^{2}} \left( V_{i+1,j} - 2 V_{i,j} + V_{i-1,j} + V_{i,j+1} - 2 V_{i,j} + V_{i,j-1} \right) = - \rho_{i,j} \quad i = 1, ..., N; ~ j = 1, ... , M
\end{equation}
In this form, given a desity $\rho$, and appropriate boundary conditions (e.g., values of $V$ on the edges of the rectangular region implict in the discretization given), solving the Poisson equation is recast as a linear problem: $Ax=b$. $A$ is a sparse matrix, $x$ is a vector of length $M \times N$, and $b$ contains the given density $\rho$.
Write a code that would use an iterative method of your choice to solve Poisson's equation. Imagine that the charge density is given, and that the boundary consitions are also specified.
!!!!2. Numerical solution of Shrodinger's equation
Consider a one dimensional potential well that is a truncated harmonic oscillator potential:
\begin{equation}
V(x) = \frac{1}{2} x^{2} \quad 0 < x < 1, \quad\quad V = \infty \quad x=0,1
\end{equation}
Take as basis functions,
\begin{equation}
u_{i}(x) = \sin (q_{i} x) \quad q_{i} = i \pi
\end{equation}
Consider the Hamiltonian (in arbitrary units),
\begin{equation}
H = - \frac{d^{2} }{dx^{2}} + V(x)
\end{equation}
We wish to numerically find the energy levels of this system. First, compute the necessary matrix elements:
\begin{equation}
H_{ji} = \langle u_{j} | H | u_{i} \rangle = \int_{0}^{1} u_{j}(x) H u_{i}(x) dx
\end{equation}
Given the normalization,
\begin{equation}
S_{ji} = \langle u_{j} | u_{i} \rangle = \int_{0}^{1} u_{j}(x) u_{i}(x) dx
\end{equation}
then, Schrodinger's equation turns into an eigenvalue problem of the form:
\begin{equation}
H \alpha = E S \alpha
\end{equation}
where $\alpha$ are the components of $\psi$ in the basis given $\psi = \sum_{i} \alpha_{i} u_{i}(x)$.
Use a Lapack routine to find the energy levels of this potential. Make a judicious choice about the number of eigenfunctions you use in the diagonalization.
{{justifyright{
''Due Friday, October 7, 2016''
}}}
!!! 1. The Sucessive Over Relaxation method (SOR)
# Show that the ~Gauss-Seidel iterative method seen in class can be written as $$ x_{i}^{(k+1)} = \frac{- \sum_{j=1}^{i-1} A_{ij} x_{j}^{(k+1)}-\sum_{j=i+1}^{n} A_{ij} x_{j}^{(k)} + b_{i}}{A_{ii}}, \quad i=1, \ldots, n $$
where the residual to be made as small as possible has been defined as $r_{i}^{(k)} = x_{i}^{(k+1)} - x_{i}^{(k)}$.
# "Over relaxation" is introduced in order to speed up convergence of the iteration, and consists in introducing a new iteration that overcorrects the residual,$$ x_{i}^{(k+1)} = x_{i}^{(k)} + \omega r_{i}^{(k)}$$ where $\omega$ is known as the relaxation parameter. Show that if we write, in matrix form, $x^{(k+1)} = B x^{(k)} + C$, then for this new method we have for the iteration $$ B_{\rm SOR} = \left( \mathbb{1} + \omega D^{-1}L \right)^{-1} \left[ (1-\omega)\mathbb{1} - \omega D^{-1}U \right]$$
with the same notation used in class for $D, L$ and $U$.
# Show that the determinant of $B$ is $\det (B) = (1-\omega)^{n}$.
# To estimate the spectral radius (the region of convergence of the method), assume that $B$ can be diagonalized, with eigenvalues $\lambda_{1}, \lambda_{2}, \ldots \lambda_{n}$. With this representation show that $\| B \| \ge |1-\omega|$. What is the range of $\omega$ over which this procedure will converge ?
!!! 2. Solution of a system of equations by iterative methods
Consider the system of equations $A \mathbf{x} = \mathbf{b}$ where $$ A = \left( \begin{array}{cccc} 4 & -1 & -1 & 0 \\ -1 & 4 & 0 & -1 \\ -1 & 0 & 4 & -1 \\ 0 & -1 & -1 & 4 \end{array} \right), \quad\quad b = \left( \begin{array}{c} 1 \\ 2 \\ 0 \\1 \end{array} \right). $$
Write a code that solves this system by using (i) the Jacobi method, (ii) the ~Gauss-Seidel method, and (iii) the SOR method. From the output of the iterations, compare empirically the rate of convergence of the three methods.
!!!3. Solution of Poisson equation in two dimensions
Poisson's equation relates the electrostatic potential $V(x,y)$ to a charge distribution $\rho(x,y)$. In arbitrary units, and in two spatial dimensions we write,
\begin{equation}
\frac{\partial^{2} V(x,y)}{\partial x^{2}} + \frac{\partial^{2} V(x,y)}{\partial y^{2}} = - \rho(x,y)
\end{equation}
By some specific discretization scheme, the equation is written as (with $V_{ij} = V(x_{i},y_{j})$, $\Delta x = \Delta y$):
\begin{equation}
\frac{1}{\Delta x^{2}} \left( V_{i+1,j} - 2 V_{i,j} + V_{i-1,j} + V_{i,j+1} - 2 V_{i,j} + V_{i,j-1} \right) = - \rho_{i,j} \quad i = 1, ..., N; ~ j = 1, ... , M
\end{equation}
In this form, given a density $\rho$, and appropriate boundary conditions (e.g., values of $V$ on the edges of the rectangular region implicit in the discretization given), you can think of solving the Poisson equation as solving a linear problem: $Ax=b$. $A$ is a sparse matrix, $x$ is a vector of length $M \times N$, and $b$ contains the given density $\rho$. Alternatively, you may also directly write an iterative procedure to solve the discrete system of $N \times M$ equations for the potential $V_{ij}$.
Write a code that would use an iterative method of your choice to solve Poisson's equation. Imagine that the charge density is given, and that the boundary conditions are also specified.
{{justifyright{
''Due Tuesday, October 2, 2018''
}}}
!!! Question 1. Spline interpolation with explicit solution of a tridiagonal system
Spline polynomials are smooth, piecewise polynomials that are extensively used in Physics to interpolate a sample of discrete points.
We have found experimentally the pressure $p$ as a function of the molar volume $v$ of a binary mixture at constant temperature and composition (in some arbitrary units). The results are shown in the following table.
|v | 1.0 | 1.1 | 1.2 | 1.3 | 1.4 | 1.5 |
| p | 4.7 | 3.5 | 3.0 | 2.7 | 3.0 | 2.4 |
# Write your own code to find a cubic spline interpolation polynomial $p=p(v)$ of this data in the interval given. Plot the data along with the spline interpolation on the same graph.
# Use the spline interpolation to find the value of the volume at which $p=3.25$. Do it approximately on the graph first.
# Compute the value of the volume at which $p=3.25$ by solving the equation $f(v)=p(v)-3.25=0$ by using the ~Newton-Raphson method ''WITH'' the interpolated spline function.
# Additional information
## Documentation about spline polynomials can be found at Wolfram [[MathWorld|http://mathworld.wolfram.com/CubicSpline.html]]. Eq. (1) defines the collection of polynomials, and Eqs. (6-9) define the mapping between these coefficients and the set $D_{i}$. The latter is the one that is actually calculated in practice by solving the tridiagonal system in Eq. (18).
## Use the algorithm given in class to solve tridiagonal equations.
!!! Question 2. Spline interpolation. Routines provided
Rewrite the code of question 1 but now used spline interpolation routines used by others.
# If you use FORTRAN, you can use the routines spline.f to compute the spline coefficients and seval.f to evaluate the polynomials at interpolated points. The routines are provided in the class github repository.
# If you use C++, you can use the functions spline.c to compute the polynomial coefficients, splint.c to evaluate the polynomials obtained. You will also need to add in the compilation nrutil.c and the header file nrutil.h. All of them are provided in the class github repository.
!!! Question 3. Numerical solution of Schroedinger's equation
Consider a one dimensional potential well that is a truncated harmonic oscillator potential:
\begin{equation}
V(x) = \frac{1}{2} x^{2} \quad 0 < x < 1, \quad\quad V = \infty \quad x=0,1
\end{equation}
Take as basis functions,
\begin{equation}
u_{i}(x) = \frac{1}{\sqrt{2}} \sin (q_{i} x) \quad q_{i} = i \pi
\end{equation}
Consider the Hamiltonian (in arbitrary units),
\begin{equation}
H = - \frac{d^{2} }{dx^{2}} + V(x)
\end{equation}
We wish to numerically find the energy levels of this system. First, compute the necessary matrix elements by using the trapezoidal rule,
\begin{equation}
H_{ji} = \langle u_{j} | H | u_{i} \rangle = \int_{0}^{1} u_{j}(x) H u_{i}(x) dx .
\end{equation}
Schroedinger's equation turns into an eigenvalue problem of the form:
\begin{equation}
H \alpha = E \alpha
\end{equation}
where $\alpha$ are the components of $\psi$ in the basis given $\psi = \sum_{i} \alpha_{i} u_{i}(x)$.
Use a Lapack routine to find the energy levels of this potential. Make a judicious choice about the number of eigenfunctions you use in the diagonalization.
{{justifyright{
''Due Wednesday, October 7, 2020''
}}}
!!! 0. Format requirements
In order to streamline uploading into Canvas and grading, please adhere to the following:
# Submit one separate source file for each of the three questions. If there is a main() and functions for a question, please include them in the same file.
# If you submit a Jupyter notebook, please also submit a .pdf export to simplify quick referencing.
# Results, plots, and discussion should be contained in a separate, single file (any format that you wish. If the format is obscure, please export to .pdf and submit the pdf file).
# Do not submit executables (after compilation). They are not portable in principle, and cannot be used.
!!! 1. Direct and inefficient direct matrix inversion to solve a differential equation
We will use the LAPACK routine {{{dgesv}}} to solve a partial differential equation.
The heat diffusion equations in one dimension is,
$$ \frac{\partial u(x,t)}{\partial t} = D \frac{\partial^{2} u}{\partial x^{2}}$$
Assume we denote by $u_{i}^{n}=u(x_{i},t_{n})$ the solution at a discrete set of points in some interval of the $x$ axis, $x_{i} = i \Delta x$. $i = 0, \ldots N$, and corresponding to some discrete time $t_{n} = n \Delta t$, with $\Delta x$ and $\Delta t$ the discrete increments in space and time respectively. We are told that the differential equation can be written in these discretized variables as (we will study this later in the semester)
$$ \frac{1}{\Delta t} \left( u_{i}^{n+1}-u_{i}^{n} \right) = \frac{D}{(\Delta x)^{2}} \left( u_{i+1}^{n+1} - 2 u_{i}^{n+1} + u_{i-1}^{n+1} \right) $$
By rearranging the terms of this equation, an iterative solution can be constructed that relates the solution at some time $t_{n}$, $u_{i}^{n}$, to the solution at the next time step $t_{n+1}$, $u_{i}^{n+1}$. The iteration is,
$$ -a u_{i-1}^{n+1} + (1+2a) u_{i}^{n+1} - a u_{i+1}^{n+1} = u_{i}^{n}$$
with $a = D \Delta t/(\Delta x)^{2}$. This iteration can be written in matrix form $Ax=b$ where $x=\{u_{i}^{n+1} \}$, $b=\{u_{i}^{n} \}$ and the matrix $A$ is $$A = \left( \begin{array}{cccccc} 1+2a & -a & 0 & \ldots& & 0 \\ -a & 1+2a & -a & 0 & \dots & 0 \\ 0 & -a & 1+2a & -a & \ldots & 0 \\ & \ldots & \ldots & \ldots & & \\ 0 & \ldots & 0 & & -a & 1+2a \end{array} \right)$$
As boundary conditions at the end points, we take fixed $u$ for all times: $u_{0} = 0 $ and $u_{N} = 0$.
Write a code that uses the LAPACK routine {{{dgesv}}} to solve the linear system of equations given which amounts to solving the partial differential equation iteratively: Solve for a number of time steps $n = 1, 2, \ldots$ with an initial condition at $t=0$ $u_{k}^{0}=1$ for $k = N/2$ and the remaining $u_{k}^{0} = 0$. In Python, you will use {{{dgesv}}} indirectly through the use of {{{numpy.linalg.solve}}}.
Plot the solution for a few time steps as it evolves in time.
!!! 2. Numerical solution of the stationary Shrodinger's equation
Consider a one dimensional potential well that is a truncated harmonic oscillator potential:
\begin{equation}
V(x) = \frac{1}{2} x^{2} \quad 0 < x < 1, \quad\quad V = \infty \quad x=0,1
\end{equation}
Take as basis functions,
\begin{equation}
u_{i}(x) = \sin (q_{i} x) \quad q_{i} = i \pi
\end{equation}
Consider the Hamiltonian (in arbitrary units),
\begin{equation}
H = - \frac{d^{2} }{dx^{2}} + V(x)
\end{equation}
We wish to numerically find the energy levels of this system. First, compute the necessary matrix elements:
\begin{equation}
H_{ji} = \langle u_{j} | H | u_{i} \rangle = \int_{0}^{1} u_{j}(x) H u_{i}(x) dx
\end{equation}
Given the normalization,
\begin{equation}
S_{ji} = \langle u_{j} | u_{i} \rangle = \int_{0}^{1} u_{j}(x) u_{i}(x) dx
\end{equation}
then, Schrodinger's equation turns into an eigenvalue problem of the form:
\begin{equation}
H \alpha = E S \alpha
\end{equation}
where $\alpha$ are the components of $\psi$ in the basis given $\psi = \sum_{i} \alpha_{i} u_{i}(x)$.
Use a Lapack routine to find the energy levels of this potential. Make a judicious choice about the number of eigenfunctions you use in the diagonalization. In Python you may use {{{numpy.linalg.eig}}}.
!!! 3. Convergence of numerical differentiation
Write code that verifies the relatively faster convergence of a numerical derivative computed by central differences compared to a derivative computed by forward differences. Evaluate the lowest order forward and central difference approximation to the derivative of $f(x) = \tan (x)$ in the interval $[-2,2]$ for $n$ equally spaced points. Compare the results to the exact values of the derivative. For a given $n$ find the maximum error on the interval from each method, defining the error as the absolute value of the ratio of the difference between numerical and analytic derivatives. Then make a log-log plot of this error versus $n$ over at least one order of magnitude in $n$ (say from $n = 100$ to $n=5000$) to show that the forward difference estimate converges linearly with $n$, while the central difference estimate converges as $n^{2}$.
{{justifyright{
''Due Tuesday, October 14, 2014''
}}}
# Write code that verifies the relatively faster convergence of a numerical derivative computed by central differences compared to a derivative computed by forward differences. Evaluate the lowest order forward and central difference approximation to the derivative of $f(x) = \tan (x)$ in the interval $[-2,2]$ for $n$ equally spaced points. Compare the results to the exact values of the derivative. For a given $n$ find the maximum error on the interval from each method, defining the error as the absolute value of the ratio of the difference between numerical and analytic derivatives. Then make a log-log plot of this error versus $n$ over at least one order of magnitude in $n$ (say from $n = 100$ to $n=5000$) to show that the forward difference estimate converges linearly with $n$, while the central difference estimate converges as $n^{2}$.
# Numerical intergation over an unbounded domain. The integration methods discussed in class consider definite integrals over a finite domain $(a,b)$. However, it is possible to either transform the integrands or to analyze convergence with the limits of integration to evaluate integrals over unbounded domains. Propose a suitable plan, and write the associate code to evaluate the following integrals for $s = 0.5, 0.6, 0.7, ..., 3.0$.
## $$\int_{0}^{\infty} \left( x^{3} + sx \right)^{-1/2} dx,$$
## $$\int_{0}^{\infty} \left( x^{2} + 1 \right) ^{-1/2} e^{-sx} dx$$
# Numerical integration of a predator-prey model. Consider two species of fish that coexist in the same lake. Let $x$ be the number of small fish of type A which only eat grass from the bottom, and $y$ the number of large fish of type B which only eat small fish of type A, but no grass. Alone, the population of A might be expected to grow exponentially, whereas B alone (no A to eat) would die exponentially. However, the two species interact: whenever a fish B encounters a fish A, it will eat it. This interaction contributes to the increase in the population of B and a decrease in A, both rates being proportional to the quantity $xy$. A general mathematical model of this system is $$ \frac{dx}{dt} = \epsilon_{A} x - \gamma_{A} xy $$
$$ \frac{dy}{dt} = - \epsilon_{B} y + \gamma_{B} xy $$ where $\epsilon_{A}$ is the growth rate of A alone, $\epsilon_{B}$ the extinction of B alone, and $\gamma_{A}$ and $\gamma_{B}$ the interaction coefficients. All four are assumed to be constant. Write a code that integrates this system of differential equations with $\epsilon_{A} = \epsilon_{B} = \gamma_{A} = \gamma_{B} =1$, and initial conditions at $t=0$ $x=0.1$ and $y=0.1$.
## Integrate the system up to $t=10$.
## Plot $y$ versus $x$ and provide an interpretation of the resulting graph.
[[pdf version|http://homepages.spa.umn.edu/~vinals/tspot_files/phys4041/2014/PHYS4041_Homework4.pdf]].
{{justifyright{
''Due Friday, October 14, 2016''
}}}
!!!1. Numerical solution of Shrodinger's equation
Consider a one dimensional potential well that is a truncated harmonic oscillator potential:
\begin{equation}
V(x) = \frac{1}{2} x^{2} \quad 0 < x < 1, \quad\quad V = \infty \quad x=0,1
\end{equation}
Take as basis functions,
\begin{equation}
u_{i}(x) = \sin (q_{i} x) \quad q_{i} = i \pi
\end{equation}
Consider the Hamiltonian (in arbitrary units),
\begin{equation}
H = - \frac{d^{2} }{dx^{2}} + V(x)
\end{equation}
We wish to numerically find the energy levels of this system. First, compute the necessary matrix elements:
\begin{equation}
H_{ji} = \langle u_{j} | H | u_{i} \rangle = \int_{0}^{1} u_{j}(x) H u_{i}(x) dx
\end{equation}
Given the normalization,
\begin{equation}
S_{ji} = \langle u_{j} | u_{i} \rangle = \int_{0}^{1} u_{j}(x) u_{i}(x) dx
\end{equation}
then, Schrodinger's equation turns into an eigenvalue problem of the form:
\begin{equation}
H \alpha = E S \alpha
\end{equation}
where $\alpha$ are the components of $\psi$ in the basis given $\psi = \sum_{i} \alpha_{i} u_{i}(x)$.
Use a Lapack routine to find the energy levels of this potential. Make a judicious choice about the number of eigenfunctions you use in the diagonalization.
!!!2. Convergence of numerical differentiation
Write code that verifies the relatively faster convergence of a numerical derivative computed by central differences compared to a derivative computed by forward differences. Evaluate the lowest order forward and central difference approximation to the derivative of $f(x) = \tan (x)$ in the interval $[-2,2]$ for $n$ equally spaced points. Compare the results to the exact values of the derivative. For a given $n$ find the maximum error on the interval from each method, defining the error as the absolute value of the ratio of the difference between numerical and analytic derivatives. Then make a log-log plot of this error versus $n$ over at least one order of magnitude in $n$ (say from $n = 100$ to $n=5000$) to show that the forward difference estimate converges linearly with $n$, while the central difference estimate converges as $n^{2}$.
!!!3. Numerical integration over an unbounded domain.
The integration methods discussed in class consider definite integrals over a finite domain $(a,b)$. However, it is possible to either transform the integrands or to analyze convergence with the limits of integration to evaluate integrals over unbounded domains. Propose a suitable plan, and write the associate code to evaluate the following integrals for $s = 0.5, 0.6, 0.7, ..., 3.0$.
## $$\int_{0}^{\infty} \left( x^{3} + sx \right)^{-1/2} dx,$$
## $$\int_{0}^{\infty} \left( x^{2} + 1 \right) ^{-1/2} e^{-sx} dx$$
{{justifyright{
''Due Tuesday, October 9, 2018''
}}}
!!! 1. Solution of a system of equations by iterative methods
Consider the system of equations $A \mathbf{x} = \mathbf{b}$ where $$ A = \left( \begin{array}{cccc} 4 & -1 & -1 & 0 \\ -1 & 4 & 0 & -1 \\ -1 & 0 & 4 & -1 \\ 0 & -1 & -1 & 4 \end{array} \right), \quad\quad b = \left( \begin{array}{c} 1 \\ 2 \\ 0 \\1 \end{array} \right). $$
Write a code that solves this system by using (i) the Jacobi method, and (ii) the ~Gauss-Seidel method. From the output of the iterations, compare empirically the rate of convergence of the two methods.
!!! 2. Solution of an integral equation
We need to find the solution of the following integral equation for $f(x)$ in the interval $(0,1)$,$$
\int_{0}^{1} dx e^{-(x-y)^{2}/\epsilon^{2}} \frac{df(x)}{dx} = 2 \epsilon \sqrt{\pi} ( y - 3 y^{2} + 2 y^{3}),
$$
where the function $f(x)$ needs to satisfy the boundary conditions $f(0)=f(1)=0$, and $\epsilon = 0.01$.
Discretize $x_{j}= j h, j = 0, \ldots n$ (take $n = 101$) for some $h$, so that $x_{0} = 0$ and $x_{n} = 1$, and $y_{i} = i h, i = 0, \ldots n$. Consistent with this discretization, one defines $f_{j} = f(x_{j})$. The boundary conditions are then $f_{0} = f_{n} = 0$.
# Write the discretized version of the integral equation by using the trapezoidal rule for the integral and central differences for the derivative at the inner points, the forward difference at $x=0$ and the backward difference at $x=1$.
# Use the discretized version to write the code that will find the solution $f(x)$ in $(0,1)$.
{{justifyright{
''Wednesday, October 14, 2020''
}}}
!!! 1. Extrapolation
A physical quantity $x$ is assumed to depend on the pressure of a gas $p$ according to
$$
x = c_{0} + c_{2} p^{2} + c_{3} p^{3} + c_{6} p^{6}
$$
where $c_{0}, c_{2}, c_{3}$ and $c_{6}$ are nonzero constants. By using repeated extrapolation, determine the value of $x$ in a vacuum (zero pressure) from the following measurements,
|! p (mmHg) | 0.8 | 0.4 | 0.2 | 0.1 | 0.05 |
|! x | 740 | 487 | 475 | 485 | 489 |
!!! 2. Numerical integration of a predator-prey model
Consider two species of fish that coexist in the same lake. Let $x$ be the number of small fish of type A which only eat grass from the bottom, and $y$ the number of large fish of type B which only eat small fish of type A, but no grass. Alone, the population of A might be expected to grow exponentially, whereas B alone (no A to eat) would die exponentially. However, the two species interact: whenever a fish B encounters a fish A, it will eat it. This interaction contributes to the increase in the population of B and a decrease in A, both rates being proportional to the quantity $xy$. A general mathematical model of this system is $$ \frac{dx}{dt} = \epsilon_{A} x - \gamma_{A} xy $$
$$ \frac{dy}{dt} = - \epsilon_{B} y + \gamma_{B} xy $$ where $\epsilon_{A}$ is the growth rate of A alone, $\epsilon_{B}$ the extinction of B alone, and $\gamma_{A}$ and $\gamma_{B}$ the interaction coefficients. All four are assumed to be constant. Write a code that integrates this system of differential equations with $\epsilon_{A} = \epsilon_{B} = \gamma_{A} = \gamma_{B} =1$, and initial conditions at $t=0$ $x=0.1$ and $y=0.1$.
# Integrate the system up to $t=10$.
# Plot $y$ versus $x$ and provide an interpretation of the resulting graph.
!!! 2. The SEIR model of viral infection
Epidemiological modeling relies on the so called "compartmental models". A widely used one is the SEIR model ("Susceptible, Exposed, Infected, and Recovered"). The model is defined by the following set of [[differential equations|https://en.wikipedia.org/wiki/Compartmental_models_in_epidemiology#The_SEIR_model]],
\begin{equation}
\frac{dS}{dt} = \mu N - \mu S -\beta \frac{IS}{N}
\end{equation}
\begin{equation}
\frac{dE}{dt} = \beta \frac{IS}{N} - (\mu + a)E
\end{equation}
$$
\frac{dI}{dt} = aE - (\gamma+ \mu) I
$$
$$
\frac{dR}{dt} = \gamma I - \mu R
$$
so that $S+E+I+R = N$, constant. $S$ is the number of susceptible individuals as a function of time, $E$ the number of exposed, $I$ the number of infected, and $R$ the number of recovered. $N$ is the total population, constant in this version of the model. The various coefficients are rates (e.g., $\beta$ is the rate of infection, etc.).
This model has a basic parameter
$$
R_{0} = \frac{a}{\mu + a} \; \frac{\beta}{\mu + \gamma}.
$$
If $R_{0} < 1$ this model has a fixed point in which everyone is healthy (as $t \rightarrow \infty, S = N, E=I=R=0$). If on the other hand $R_{0}>1$ there is another fixed point as $t \rightarrow \infty$ in which infection persists.
Write a code that will solve this system of ordinary differential equations for arbitrary initial conditions, and verify the asymptotic convergence to a fixed point at long times. Plot the results.
You may use the following values of the coefficients: $\mu = 0.1, a = 0.1, \gamma = 0.1$ and $\beta = 0.39, 0.41$.
{{justifyright{
''Due Tuesday, November 4/11, 2014''
}}}
!!!!1. ~Lax-Wendroff method. Due November 4.
Write a program that uses the ~Lax-Wendroff method to solve the PDE:
$$ \frac{\partial u}{\partial t} = - v \frac{\partial u}{\partial x}$$
with $v$ a given constant. As given in class, the method as applied to this equation is
$$u_{j}^{n+1}=u_{j}^{n}-\frac{v \Delta t}{\Delta x} \left( u_{j+1/2}^{n+1/2} - u_{j-1/2}^{n+1/2} \right)$$
with
$$u_{j+1/2}^{n+1/2} = \frac{1}{2} \left( u_{j+1}^{n} + u_{j}^{n} \right) - \frac{v \Delta t}{2 \Delta x} \left( u_{j+1}^{n} - u_{j}^{n} \right) $$
and a similar definition for $u_{j-1/2}^{n+1/2}$. Space is discretized by $\Delta x$ and $j$ is the spatial index, whereas time has been discretized by $\Delta t$ with corresponding index $n$.
Let $a = |v|\Delta t/\Delta x \le 1$ denote the Courant condition for stability.
Consider as initial condition $u(x,t=0) = \sin (2 \pi x)$ for $0 \le x \le 1$. If your discretization is such that $x_{0} = 0$ and $x_{N} = 1$, have the code solve the equation with "periodic boundary conditions": $u_{N} = u_{0}$ at each time (and should you need it $u_{N+1} = u_{1}$, etc.). For the purposes of the homework, choose $N = 256$ and $v = 1$.
Run your code with $a = 0.99$. Since $v = 1$, the solution is advected in such a way that the solution should exactly lie on top of the initial condition for times $t = 1, 2, 3, ...$. Confirm this graphically for $t = 1, 2$.
Run the code until $t = 0.1$ using both $a = 0.99$ and $a = 1.1$. Plot the solutions and comment on any differences.
!!!!2. Steady state temperature field in a porous medium. Due November 11.
The temperature distribution in a porous rock in the presence of ground water flow directed along the $x$ axis is described by the equation
$$ v \frac{\partial T}{\partial x} = \kappa \nabla^{2} T$$
where $v$ is the constant water velocity and $\kappa$ is the thermal diffusivity coefficient. The domain of solution is the rectangle $0 < xv/\kappa < 30$ $0 < yv/\kappa < 10$.
Compute the temperature distribution using the Alternating Direction Implicit method combined with a finite difference method of your choice. The boundary conditions $T = 65$ at $x=0$, $T = 25$ at $x = 30 \kappa/v$, $T = 25$ at $y = 10 \kappa/v$ and $\partial T/\partial y = 0$ at $y=0$.
!!!! Notes on the ADI method
It is useful to be able to write the system of equations of Problem 2 as the usual system $Ax = b$ where both $x$ and $b$ are arrays. Consider for example the two dimensional Poisson equation
$$ T_{i+1 j} - 2 T_{ij} + T_{i11 j} + T_{i j+1} - 2 T_{ij} + T_{i j-1} = - \Delta x^{2} \rho_{ij}$$
for $i = 1, ... , N$ and $j = 1, ..., M$.
We can now construct a long vector $\mathbf{v}$ of length $N \times M$ by linking together the @@color(red):rows@@ of the matrix $T_{ij}$:
$$ v_{r} = T_{ij} \quad \quad r = (i-1)M + j$$
This mapping can be inverted as
$$ i = 1+ {\rm INT}\left( \frac{r-1}{M} \right) \quad\quad j = 1 + (r-1){\rm MOD}(M)$$
The Poisson equation then transforms to
$$ v_{r-M} + v_{r-1} - 4 v_{r} + v_{r+1} + v_{r+M} = - \Delta x^{2} \rho_{r}$$ which can now written in the standard form $Ax = b$.
For the ADI method we need to construct another long vector $\mathbf{w}$ by linking together the @@color(red):columns@@ of the matrix $T_{ij}$:
$$ w_{s} = T_{ij} \quad\quad s = (j-1)N + i $$
which can in turn be inverted as
$$ j = 1 + {\rm INT}\left( \frac{s-1}{N} \right) \quad\quad i = 1 + (s-1){\rm MOD}(N) $$
With the aid of @@color(red):both@@ vectors, the Poisson equation can be written as
$$ w_{s+1} - 2w_{s} + w_{s-1} + v_{r+1} - 2 v_{r} + v_{r-1} = - \Delta x^{2} \rho_{ij}$$
which one can write in matrix form as
$$ A_{1} \mathbf{v} + A_{2} \mathbf{w} = b$$
so that $A_{1}$ only acts on the rows of $T_{ij}$ and $A_{2}$ on the columns. We also define the reordering matrix $U$ as
$$ \mathbf{w} = U \mathbf{u}$$
In this notation, the ADI iteration can be written as
\begin{eqnarray}
\left( A_{1} + I \right) \mathbf{v}^{n+1/2} & = & b - (A_{2} \mathbf{w}^{n} - \mathbf{v}^{n} ) \\
\mathbf{w}^{n+1/2} & = & U \mathbf{v}^{n+1/2} \\
\left( A_{2} + I \right) \mathbf{w}^{n+1} & = & b - (A_{1} \mathbf{v}^{n+1/2} - \mathbf{w}^{n+1/2} )
\end{eqnarray}
where $I$ is the identity matrix.
{{justifyright{
''Due Friday, November 4, 2016''
}}}
!!! ~Lax-Wendroff method.
Write a program that uses the ~Lax-Wendroff method to solve the PDE:
$$ \frac{\partial u}{\partial t} = - v \frac{\partial u}{\partial x}$$
with $v$ a given constant. As given in class, the method as applied to this equation is
$$u_{j}^{n+1}=u_{j}^{n}-\frac{v \Delta t}{\Delta x} \left( u_{j+1/2}^{n+1/2} - u_{j-1/2}^{n+1/2} \right)$$
with
$$u_{j+1/2}^{n+1/2} = \frac{1}{2} \left( u_{j+1}^{n} + u_{j}^{n} \right) - \frac{v \Delta t}{2 \Delta x} \left( u_{j+1}^{n} - u_{j}^{n} \right) $$
and a similar definition for $u_{j-1/2}^{n+1/2}$. Space is discretized by $\Delta x$ and $j$ is the spatial index, whereas time has been discretized by $\Delta t$ with corresponding index $n$.
Let $a = |v|\Delta t/\Delta x \le 1$ denote the Courant condition for stability.
Consider as initial condition $u(x,t=0) = \sin (2 \pi x)$ for $0 \le x \le 1$. If your discretization is such that $x_{0} = 0$ and $x_{N} = 1$, have the code solve the equation with "periodic boundary conditions": $u_{N} = u_{0}$ at each time (and should you need it $u_{N+1} = u_{1}$, etc.). For the purposes of the homework, choose $N = 256$ and $v = 1$.
Run your code with $a = 0.99$. Since $v = 1$, the solution is advected in such a way that the solution should exactly lie on top of the initial condition for times $t = 1, 2, 3, ...$. Confirm this graphically for $t = 1, 2$.
Run the code until $t = 0.1$ using both $a = 0.99$ and $a = 1.1$. Plot the solutions and comment on any differences.
{{justifyright{
''Due Tuesday, October 30, 2018''
}}}
!!! 1. Eigenvalue problem for a differential equation
Write a code that computes the smallest eigenvalue $\lambda$ for which the problem $$
\frac{d}{dx} \left[ (1+x^{2}) \frac{dy}{dx} \right] + \lambda y = 0, \quad\quad y(-1)=y(1)=0
$$
has a nontrivial solution.
In order to do it, you can use one of the following two methods. Either follow the method described in class in which a new variable $z(x)=\lambda$ is introduced, with $dz/dx = 0$, or user a matrix method in the following way: discretize the equation as, $$
\frac{d}{dx} \left[ p(x) \frac{dy}{dx} \right] = \frac{1}{h} \left[ p_{n+1/2}\left( \frac{y_{n+1}-y_{n}}{h} \right) - p_{n-1/2} \left( \frac{y_{n}-y_{n-1}}{h} \right) \right] + c_{1}(x) h^{2} + c_{2} h^{4} + \ldots ,
$$
where we have also included the asymptotic error of the discretization. Once the differential equation is discretized, you can find the smallest eigenvalue of the resulting matrix.
!!! 2. ~Lax-Wendroff method
Write a program that uses the ~Lax-Wendroff method to solve the PDE:
$$ \frac{\partial u}{\partial t} = - v \frac{\partial u}{\partial x}$$
with $v$ a given constant. As given in class, the method as applied to this equation is
$$u_{j}^{n+1}=u_{j}^{n}-\frac{v \Delta t}{\Delta x} \left( u_{j+1/2}^{n+1/2} - u_{j-1/2}^{n+1/2} \right)$$
with
$$u_{j+1/2}^{n+1/2} = \frac{1}{2} \left( u_{j+1}^{n} + u_{j}^{n} \right) - \frac{v \Delta t}{2 \Delta x} \left( u_{j+1}^{n} - u_{j}^{n} \right) $$
and a similar definition for $u_{j-1/2}^{n+1/2}$. Space is discretized by $\Delta x$ and $j$ is the spatial index, whereas time has been discretized by $\Delta t$ with corresponding index $n$.
Let $a = |v|\Delta t/\Delta x \le 1$ denote the Courant condition for stability.
Consider as initial condition $u(x,t=0) = \sin (2 \pi x)$ for $0 \le x \le 1$. If your discretization is such that $x_{0} = 0$ and $x_{N} = 1$, have the code solve the equation with "periodic boundary conditions": $u_{N} = u_{0}$ at each time (and should you need it $u_{N+1} = u_{1}$, etc.). For the purposes of the homework, choose $N = 256$ and $v = 1$.
Run your code with $a = 0.99$. Since $v = 1$, the solution is advected in such a way that the solution should exactly lie on top of the initial condition for times $t = 1, 2, 3, ...$. Confirm this graphically for $t = 1, 2$.
Run the code until $t = 0.1$ using both $a = 0.99$ and $a = 1.1$. Plot the solutions and comment on any differences.
{{justifyright{
''Wednesday, October 21, 2020''
}}}
!! 0. Format requirements
In order to streamline uploading into Canvas and grading, please adhere to the following:
# Submit one separate source file for each of the three questions. If there is a main() and functions for a question, please include them in the same file.
# If you submit a Jupyter notebook, please also submit a .pdf export to simplify quick referencing.
# Results, plots, and discussion should be contained in a separate, single file (any format that you wish. If the format is obscure, please export to .pdf and submit the pdf file).
# Do not submit executables (after compilation). They are not portable in principle, and cannot be used.
!! 1. Boundary value problem: solve the Schrödinger equation by matching
We wish to solve the second order differential equation,
$$
- \frac{\hbar^{2}}{2m} \frac{d^{2} \psi(x)}{dx^{2}} + V(x) \psi(x) = E \psi(x)
$$
with
$$
V(x) = \left\{ \begin{array}{ll} -V_{0} = -83MeV & |x| \le a = 2 \ {\rm fm} \\ 0 & |x| \gt a = 2 \ {\rm fm} \end{array} \right.
$$
These values are typical for a nucleus.
We first introduce dimensionless variables by considering a typical nucleon mass so that $2m/\hbar^{2} = 0.0483 MeV^{-1} {\rm fm}^{-1}$. We also define the dimensionless wavenumber $k$ as
$$
k^{2} = - \frac{2m}{\hbar^{2}} E.
$$
//In dimensionless variables, the equations to solve are //
\begin{eqnarray}
\frac{d^{2} \psi}{dx^{2}} + \left( \frac{2m}{\hbar^{2}} V_{0} - k^{2} \right) \psi(x) &=& 0, \quad |x| \le a \\
\frac{d^{2} \psi}{dx^{2}} - k^{2} \psi(x) &= & 0, \quad |x| > a
\end{eqnarray}
//We need next to specify boundary conditions//. This is somewhat complicated here. If the particle is inside the potential well $(-a,a)$, then we know that the wave function decays exponentially for both $x \pm \infty$. We write this result for both the left and right domains as
\begin{eqnarray}
\psi(x) & \rightarrow & \psi_{R}(x) = e^{-kx}, & x & \rightarrow + \infty \\
\psi(x) &\rightarrow & \psi_{L}(x) = e^{kx}, & x & \rightarrow - \infty
\end{eqnarray}
''The numerical problem involves finding the value of $k^{2}$ (proportional to the energy $E$) that gives a continuum function smoothly connecting $\psi_{L}$ to $\psi_{R}$.''
It is suggested that you proceed as follows.
# Define a far left and far right point $|\pm x_{max}| \gg a$. Also define a matching point $x_{m}$. In principle they are all arbitrary and the solution should not depend on this choice. Choose $x_{m} \gt \approx a$. Also choose an initial guess for $k$ (so that $E > - V_{0}$).
# Use any ODE solver to step $\psi_{L}(x)$ from $-x_{max}$ towards the origin, until you reach $x_{m}$, with a value obtained $\psi_{L}(x_{m})$.
# Repeat the process starting from $+x_{max}$ to find an approximate value $\psi_{R}(x_{m})$.
# We now that the wave function and its derivative have to be continuous everywhere, including at the matching point $x_{m}$. The two solutions that you will find with some arbitrary starting values at $x = \pm x_{max}$, will not satisfy this property. We will define an error based on the ratio $\psi^{\prime}/\psi$, $$
\epsilon = \frac{\psi_{L}^{\prime}(x_{m})/\psi_{L}(x_{m}) - \psi_{R}^{\prime}(x_{m})/\psi_{R}(x_{m})}{\psi_{L}^{\prime}(x_{m})/\psi_{L}(x_{m}) + \psi_{R}^{\prime}(x_{m})/\psi_{R}(x_{m})}
$$
If both the function and the derivative are continuous, then $\epsilon = 0$.
# For your initial guess of $k$, the error $\epsilon$ will not be zero. Devise an iterative procedure for $k$ and repeat steps 2-4 until some tolerance is reached for the error $\epsilon$.
# Plot the final, converged, wavefunction.
!! 2. Parametric resonance and the Mathieu equation
Standard resonance occurs when a harmonic oscillator is driven by a periodic function of frequency close to the natural frequency of the oscillator. There is a different type of resonance when, for example, the support point of a pendulum is vibrated, or when a container with a layer of fluid is vibrated in the vertical direction (surface waves).
In some units, the variable of interest satisfies the Mathieu equation,
$$
\frac{d^{2}x}{dt^{2}} + \left( p - 2 q \cos (2 t) \right) x = 0.
$$
where $p$ and $q$ are two constants. If $q = 0$ this is the equation of the harmonic oscillator with angular frequency $\sqrt{p}$. However, in this equation the "frequency" itself is a periodic function. This simple equation cannot be solved analytically.
Use the ~Runge-Kutta method of order 4 to find the solution to this equation, under the following conditions,
# Use $p = a_{0}(q)$ and $q=1$ with initial conditions $x(0) = 1$, $x^{\prime}(0) = 0$. Discuss the result.
# Use $p = b_{1}(q)$ and $q=1$ with initial conditions $x(0) = 0$, $x^{\prime}(0) = 1$.
For the evaluation of $p$ use the results,
\begin{eqnarray}
p = a_{0}(q) & = & -\frac{1}{2} q^{2} + \frac{7}{128} q^{4} - \frac{29}{2304} q^{6} + \ldots \\
p = b_{1}(q) & = & 1 - q - \frac{1}{8}q^{2} + \frac{1}{64} q^{3} - \frac{1}{1536}q^{4} + \ldots
\end{eqnarray}
{{justifyright{
''Due Friday, November 11, 2016''
}}}
!!! Steady state temperature field in a porous medium.
The temperature distribution in a porous rock in the presence of ground water flow directed along the $x$ axis is described by the equation
$$ v \frac{\partial T}{\partial x} = \kappa \nabla^{2} T$$
where $v$ is the constant water velocity and $\kappa$ is the thermal diffusivity coefficient. The domain of solution is the rectangle $0 < xv/\kappa < 30$ $0 < yv/\kappa < 10$.
In order to find the solution, we will solve the time dependent convection-diffusion equation instead $$
\frac{\partial T}{\partial t} = - v \frac{\partial T}{\partial x} + \kappa \left( \frac{\partial^{2} T}{\partial x^{2}} + \frac{\partial^{2} T}{\partial y^{2}} \right),
$$
and iterate it in time until it reaches a steady state $\partial T/\partial t = 0$, at which point we have solved the required equation.
By using forward differences for the time derivative, central differences for the spatial derivatives, and the Alternatig Direction Implicit method to consider a semi-implicit solution in the two spatial directions, obtain the temperature distribution given the following boundary conditions are $T = 65 \;$ at $\; x=0$, $T = 25 \;$ at $\; x = 30 \kappa/v$, $T = 25 \;$ at $\; y = 10 \kappa/v \;$ and $\; \partial T/\partial y = 0 \;$ at $\; y=0$. Note, you can use $\kappa = v = 1$.
{{justifyright{
''Due Tuesday, November 6, 2018''
}}}
!!! Steady state temperature field in a porous medium.
The temperature distribution in a porous rock in the presence of ground water flow directed along the $x$ axis is described by the equation
$$ v \frac{\partial T}{\partial x} = \kappa \nabla^{2} T$$
where $v$ is the constant water velocity and $\kappa$ is the thermal diffusivity coefficient. The domain of solution is the rectangle $0 < xv/\kappa < 30$ $0 < yv/\kappa < 10$.
In order to find the solution, we will solve the time dependent convection-diffusion equation instead $$
\frac{\partial T}{\partial t} = - v \frac{\partial T}{\partial x} + \kappa \left( \frac{\partial^{2} T}{\partial x^{2}} + \frac{\partial^{2} T}{\partial y^{2}} \right),
$$
and iterate it in time until it reaches a steady state $\partial T/\partial t = 0$, at which point we have solved the required equation.
By using forward differences for the time derivative, central differences for the spatial derivatives, and the Alternatig Direction Implicit method to consider a semi-implicit solution in the two spatial directions, obtain the temperature distribution given the following boundary conditions are $T = 65 \;$ at $\; x=0$, $T = 25 \;$ at $\; x = 30 \kappa/v$, $T = 25 \;$ at $\; y = 10 \kappa/v \;$ and $\; \partial T/\partial y = 0 \;$ at $\; y=0$. Note, you can use $\kappa = v = 1$.
!!!! Notes on the ADI method
As discussed in class, it is useful to be able to write the system of equations as the usual system $Ax = b$ where both $x$ and $b$ are arrays. We discretize the two dimensional domains as follows $i = 1, ... , N$ and $j = 1, ..., M$. In the first half step (in the same sequence of two steps as discussed in class) we construct a long vector $\mathbf{v}$ of length $N \times M$ by linking together the @@color(red):columns@@ of the matrix $T_{ij}$:
$$ v_{l} = T_{ij} \quad \quad l = (j-1)N + i$$
This mapping can be inverted as
$$ i = 1+ {\rm INT}\left( \frac{l-1}{N} \right) \quad\quad j = 1 + (l-1){\rm MOD}(N)$$
For the second half step, we need to construct another long vector $\mathbf{w}$ by linking together the @@color(red):rows@@ of the matrix $T_{ij}$:
$$ w_{m} = T_{ij} \quad\quad m = (i-1)M + j$$
which can in turn be inverted as
$$ j = 1 + {\rm INT}\left( \frac{m-1}{M} \right) \quad\quad i = 1 + (m-1){\rm MOD}(M) $$.
{{justifyright{
''Due Tuesday, November 25, 2014''
}}}
!!Molecular Dynamics simulation
Write a Molecular Dynamics simulation code with the following features:
* Consider a fluid system in a ''two dimensional'' square box of size $L \times L$.
* Particles interact via a standard ~Lennard-Jones potential $$V(r) = 4 \epsilon \left( \left( \frac{\sigma}{r} \right)^{12} - \left( \frac{\sigma}{r}\right)^{6} \right)$$
where $\epsilon$ and $\sigma$ are constants and assumed known.
* Use the velocity-Verlet algorithm for time integration.
* The system will be held at constant temperature $T$. In order to accomplish that, special boundary conditions will be assumed at the walls of the box (a thermostat at constant temperature). When a particle collides with any of the walls, the tangential component of the velocity of the particle is conserved, and the incident normal velocity $v_{n}$ is updated into an outgoing normal velocity $v_{n}^{\prime}$ given by [ [[van Beijeren|http://arxiv.org/abs/1411.2983]], 2014]: $$v_{n}^{\prime} = \sqrt{ - \frac{2}{\beta m} \ln \left( 1 - {\rm exp} \left( - \frac{\beta m v_{n}^{2}}{2} \right) \right) }$$ Or course, $\beta = 1 /k_{B}T$, and $m$ is the mass of the particles, all assumed identical.
* Make sure you use appropriate dimensionless variables.
Run the code and compute the internal energy, pressure and pair correlation function at some reasonable temperature.
!!!Notes on dimensionless units
Choose $\sigma$ as the unit of length: $r^{\prime}=r/\sigma$ is the dimensionless distance, and $L^{\prime}=L/\sigma$ the dimensionless lateral dimension of the box (e.g. 10). If we introduce a dimensionless potential energy $V^{\prime}=V/\epsilon$ and a dimensionless time $t^{\prime}=t/\tau$, then Newton's Law transforms to
$$ \frac{m \sigma^{2}}{\tau^{2} \epsilon} \frac{d^{2} r^{\prime}}{dt^{\prime \; 2}} = - \frac{dV^{\prime}}{dr^{\prime}}$$
To complete the choice of dimensionless variables, we define as unit or time $\tau =\sqrt{m \sigma^{2}/\epsilon}$. This definition leads to a dimensionless form of the equation of motion.
In the boundary condition, we need the ratio $mv^{2}/k_{B}T$. From the dimensionless units of space and time, a dimensionless velocity follows $v = (\sigma/\tau) v^{\prime}$. Therefore we find
$$ \frac{mv^{2}}{k_{B}T} = \frac{m (\sigma/\tau)^{2} v^{\prime \; 2}}{k_{B}T} = \frac{\epsilon}{k_{B} T} v^{\prime}$$ If we further define the temperature scale in units of $k_{B}T$: $\epsilon^{\prime} = \epsilon/(k_{B}T)$, then we have a complete set of dimensionless equations.
I would choose $L^{\prime} = 10$, $\epsilon^{\prime} = 1$ and about $N = 100$ particles.
{{justifyright{
''Due Friday, November 18, 2016''
}}}
!!! 1. Wave equation in a string of nonuniform density and tension
Let $y(x,t)$ be the vertical displacement of an elastic string from the horizontal, as a function of the lateral coordinate $x$, and time $t$. If the density of the string is a function of position $\rho = \rho(x)$, and the local tension on the string is also a function of position $T=T(x)$, then the linear (small displacements) equation for the string displacement is
$$
\frac{d T(x)}{dx} \frac{\partial y(x,t)}{\partial x} + T(x) \frac{\partial^{2} y(x,t)}{\partial x^{2}} = \rho(x) \frac{\partial^{2} y(x,t)}{\partial t^{2}},
$$
where both $\rho(x) = \rho_{0}e^{\alpha x}$ and $T(x)=T_{0}e^{\alpha x}$ are assumed given.
Write a code that will integrate this equation in $x \in (0,1)$ with $y(0,t)=y(1,t) = 0$. Choose some initial condition that will yield interesting behavior, and plot the resulting evolution.
!!! 2. Wave packet in a quantum well
The time dependent Schroedinger equation for the complex wave function $\psi(x,t)$ is ($\hbar^{2}/2m = 1$),
$$
i \frac{\partial \psi(x,t)}{\partial t} = - \frac{\partial^{2} \psi(x,t)}{\partial x^{2}} + V(x) \psi(x,t)
$$
Consider an electron in a well which is initially localized at $x=5$ moving with momentum $k_{0} $,
$$
\Large{\psi(x,t=0) = e^{-\frac{1}{2} \left( \frac{x-5}{\sigma_{0}} \right)^{2}} e^{i k_{0} x}}
$$
# Write a code that solves for the motion of the electron in a square well with$V=0$ in $x\in (0,15)$, with boundary conditions $\psi(0,t) = \psi(15,t) = 0$. Take $\sigma_{0} = 0.5$, $\Delta x = 0.02$, $k_{0} = 17 \pi$, and $\Delta t=\Delta x^{2}/2$. Remember that the equation is complex, so you will have to start by writing its real and imaginary parts first, and then integrate the coupled equations that result.
# Is $\psi(x,t=0)$ an eigenstate of the Hamiltonian?
# Plot the probability distribution that represents the motion of the electron as a function of time.
{{justifyright{
''Due Friday, November 13, 2018''
}}}
!!! 1. Wave equation in a string of nonuniform density and tension
Let $y(x,t)$ be the vertical displacement of an elastic string from the horizontal, as a function of the lateral coordinate $x$, and time $t$. If the density of the string is a function of position $\rho = \rho(x)$, and the local tension on the string is also a function of position $T=T(x)$, then the linear (small displacements) equation for the string displacement is
$$
\frac{d T(x)}{dx} \frac{\partial y(x,t)}{\partial x} + T(x) \frac{\partial^{2} y(x,t)}{\partial x^{2}} = \rho(x) \frac{\partial^{2} y(x,t)}{\partial t^{2}},
$$
where both $\rho(x) = \rho_{0}e^{\alpha x}$ and $T(x)=T_{0}e^{\alpha x}$ are assumed given.
Write a code that will integrate this equation in $x \in (0,1)$ with $y(0,t)=y(1,t) = 0$. Choose some initial condition that will yield interesting behavior, and plot the resulting evolution.
!!! 2. Wave packet in a quantum well
The time dependent Schroedinger equation for the complex wave function $\psi(x,t)$ is ($\hbar^{2}/2m = 1$),
$$
i \frac{\partial \psi(x,t)}{\partial t} = - \frac{\partial^{2} \psi(x,t)}{\partial x^{2}} + V(x) \psi(x,t)
$$
Consider an electron in a well which is initially localized at $x=5$ moving with momentum $k_{0} $,
$$
\Large{\psi(x,t=0) = e^{-\frac{1}{2} \left( \frac{x-5}{\sigma_{0}} \right)^{2}} e^{i k_{0} x}}
$$
# Write a code that solves for the motion of the electron in a square well with$V=0$ in $x\in (0,15)$, with boundary conditions $\psi(0,t) = \psi(15,t) = 0$. Take $\sigma_{0} = 0.5$, $\Delta x = 0.02$, $k_{0} = 17 \pi$, and $\Delta t=\Delta x^{2}/2$. Remember that the equation is complex, so you will have to start by writing its real and imaginary parts first, and then integrate the coupled equations that result.
# Is $\psi(x,t=0)$ an eigenstate of the Hamiltonian?
# Plot the probability distribution that represents the motion of the electron as a function of time.
{{justifyright{
''Due Tuesday, December 9, 2014''
}}}
!!Monte Carlo simulation
Develop a Monte Carlo algorithm to simulate the equilibrium distribution of the Ising model in two dimensions. Consider a square lattice with $N \times N$ sites, in which each lattice site $(i,j)$ is occupied by a spin variable adopting the values $S{ij} = \pm 1$.
* The energy of a given configuration $\{ S_{ij} \}$ is given by, $$ E(\{ S_{ij} \}) = -J \sum_{<n.n.>} S_{ij}S_{kl}$$
where the constant $J$ is the interaction constant, and the sum extends only over nearest neighbors in the square lattice, without double counting interection bonds (there are four nearest neighbors in the square lattice: north, south, east and west of the spin in question).
* Use dimensionless variables: $J^{\prime} = J/k_{B}T$.
* Use periodic boundary conditions in both spatial directions. The following modular arithmetic may be useful. If $i$ is a given index (in either direction), then $ip1=i+1$ can be expressed as:$$ip1 = 1 + MOD(i,N)$$.
* Starting from a random initial condition, run two simulations for $J^{\prime} = 0.4, 0.5$. Use the Metropolis algorithm as discussed in class.
* Compute the average energy per spin $\langle E \rangle/N^{2}$ and the total magnetization per spin $\langle M \rangle = (1/N^{2}) \langle \sum_{i,j} S_{i,j} \rangle$.
{{justifyright{
''Due Friday, December 2, 2016''
}}}
!!Molecular Dynamics simulation
Write a Molecular Dynamics simulation code with the following features:
* Consider a fluid system in a ''two dimensional'' square box of size $L \times L$.
* Particles interact via a standard ~Lennard-Jones potential $$V(r) = 4 \epsilon \left( \left( \frac{\sigma}{r} \right)^{12} - \left( \frac{\sigma}{r}\right)^{6} \right)$$
where $\epsilon$ and $\sigma$ are constants and assumed known.
* Use the velocity-Verlet algorithm for time integration.
* The system will be held at constant temperature $T$. In order to accomplish that, special boundary conditions will be assumed at the walls of the box (a thermostat at constant temperature). When a particle collides with any of the walls, the tangential component of the velocity of the particle is conserved, and the incident normal velocity $v_{n}$ is updated into an outgoing normal velocity $v_{n}^{\prime}$ given by [ [[van Beijeren|http://arxiv.org/abs/1411.2983]], 2014]: $$v_{n}^{\prime} = \sqrt{ - \frac{2}{\beta m} \ln \left( 1 - {\rm exp} \left( - \frac{\beta m v_{n}^{2}}{2} \right) \right) }$$ Or course, $\beta = 1 /k_{B}T$, and $m$ is the mass of the particles, all assumed identical.
* Make sure you use appropriate dimensionless variables.
Run the code and compute the internal energy, pressure and pair correlation function at some reasonable temperature.
!!!Notes on dimensionless units
Choose $\sigma$ as the unit of length: $r^{\prime}=r/\sigma$ is the dimensionless distance, and $L^{\prime}=L/\sigma$ the dimensionless lateral dimension of the box (e.g. 10). If we introduce a dimensionless potential energy $V^{\prime}=V/\epsilon$ and a dimensionless time $t^{\prime}=t/\tau$, then Newton's Law transforms to
$$ \frac{m \sigma^{2}}{\tau^{2} \epsilon} \frac{d^{2} r^{\prime}}{dt^{\prime \; 2}} = - \frac{dV^{\prime}}{dr^{\prime}}$$
To complete the choice of dimensionless variables, we define as unit or time $\tau =\sqrt{m \sigma^{2}/\epsilon}$. This definition leads to a dimensionless form of the equation of motion.
In the boundary condition, we need the ratio $mv^{2}/k_{B}T$. From the dimensionless units of space and time, a dimensionless velocity follows $v = (\sigma/\tau) v^{\prime}$. Therefore we find
$$ \frac{mv^{2}}{k_{B}T} = \frac{m (\sigma/\tau)^{2} v^{\prime \; 2}}{k_{B}T} = \frac{\epsilon}{k_{B} T} v^{\prime}$$ If we further define the temperature scale in units of $k_{B}T$: $\epsilon^{\prime} = \epsilon/(k_{B}T)$, then we have a complete set of dimensionless equations.
I would choose $L^{\prime} = 10$, $\epsilon^{\prime} = 1$ and about $N = 100$ particles.
{{justifyright{
''Due Tuesday, November 20, 2018''
}}}
!!Molecular Dynamics simulation
Write a Molecular Dynamics simulation code with the following features:
* Consider a fluid system in a ''two dimensional'' square box of size $L \times L$.
* Particles interact via a standard ~Lennard-Jones potential $$V(r) = 4 \epsilon \left( \left( \frac{\sigma}{r} \right)^{12} - \left( \frac{\sigma}{r}\right)^{6} \right)$$
where $\epsilon$ and $\sigma$ are constants and assumed known.
* Use the velocity-Verlet algorithm for time integration.
* The system will be held at constant temperature $T$. In order to accomplish that, special boundary conditions will be assumed at the walls of the box (a thermostat at constant temperature). When a particle collides with any of the walls, the tangential component of the velocity of the particle is conserved, and the incident normal velocity $v_{n}$ is updated into an outgoing normal velocity $v_{n}^{\prime}$ given by [[[van Beijeren|http://arxiv.org/abs/1411.2983]], 2014]: $$v_{n}^{\prime} = \sqrt{ - \frac{2}{\beta m} \ln \left( 1 - {\rm exp} \left( - \frac{\beta m v_{n}^{2}}{2} \right) \right) }$$ Or course, $\beta = 1 /k_{B}T$, and $m$ is the mass of the particles, all assumed identical.
* Make sure you use appropriate dimensionless variables.
Run the code and compute the internal energy, the average interparticle separation, and the average velocity squared as a function of time.
!!!Note on dimensionless variables
Dimensionless variables are denoted by ', to distinguish them from the corresponding dimensional variable. Choose $\sigma$ as the unit of length: $r^{\prime}=r/\sigma$ is the dimensionless distance, and $L^{\prime}=L/\sigma$ the dimensionless lateral dimension of the box (e.g. 10). If we introduce a dimensionless potential energy $V^{\prime}=V/\epsilon$ and a dimensionless time $t^{\prime}=t/\tau$, then Newton's Law transforms to
$$ \frac{m \sigma^{2}}{\tau^{2} \epsilon} \frac{d^{2} r^{\prime}}{dt^{\prime \; 2}} = - \frac{dV^{\prime}}{dr^{\prime}}$$
To complete the choice of dimensionless variables, we define as unit or time $\tau =\sqrt{m \sigma^{2}/\epsilon}$. This definition leads to a dimensionless form of the equation of motion.
In the boundary condition, we need the ratio $mv^{2}/k_{B}T$. From the units of space and time, a dimensionless velocity follows $v = (\sigma/\tau) v^{\prime}$. Therefore we find
$$ \frac{mv^{2}}{k_{B}T} = \frac{m (\sigma/\tau)^{2} v^{\prime \; 2}}{k_{B}T} = \frac{\epsilon}{k_{B} T} v^{\prime \; 2}$$ If we further define the temperature scale in units of $k_{B}T$: $\epsilon^{\prime} = \epsilon/(k_{B}T)$, then we have a complete set of dimensionless equations.
I would choose $L^{\prime} = 10$, $\epsilon^{\prime} = 1$ and about $N = 100$ particles.
{{justifyright{
''Due Friday, December 9, 2016''
}}}
''Note'': You are not required to submit this homework. If you do by the due date given above, it will be graded and given back to you. The final average homework grade for this class will be computed from your best eight grades.
!!Monte Carlo simulation
Develop a Monte Carlo algorithm to simulate the equilibrium distribution of the Ising model in two dimensions. Consider a square lattice with $N \times N$ sites, in which each lattice site $(i,j)$ is occupied by a spin variable adopting the values $S{ij} = \pm 1$.
* The energy of a given configuration $\{ S_{ij} \}$ is given by, $$ E(\{ S_{ij} \}) = -J \sum_{<n.n.>} S_{ij}S_{kl}$$
where the constant $J$ is the interaction constant, and the sum extends only over nearest neighbors in the square lattice, without double counting interection bonds (there are four nearest neighbors in the square lattice: north, south, east and west of the spin in question).
* Use dimensionless variables: $J^{\prime} = J/k_{B}T$.
* Use periodic boundary conditions in both spatial directions. The following modular arithmetic may be useful. If $i$ is a given index (in either direction), then the next index, $ip1$ (i plus one), and the prior index $im1$ (i minus 1) can be expressed as:$$ip1 = 1 + MOD(i,N) \quad\quad im1 = 1+MOD(i+N-2,N)$$, where $MOD$ is the MOD function (MOD(a,b) in FORTRAN, a%b in C).
* Starting from a random initial condition, run two simulations for $J^{\prime} = 0.4, 0.5$. Use the Metropolis algorithm as discussed in class.
* Compute the average energy per spin $\langle E \rangle/N^{2}$ and the total magnetization per spin $\langle M \rangle = (1/N^{2}) \langle \sum_{i,j} S_{i,j} \rangle$.
{{justifyright{
''Due Tuesday, December 4, 2018''
}}}
!!! Method of coincidences to compute the entropy of a two level system
Consider a two level system with a ground state of energy 0, and an excited state of energy $\epsilon$. We have $N$ noninteracting particles. For some energy $m\epsilon$ ($m$ particles in the excited state and $N-m$ in the ground state), Statistical Mechanics tells us that the number of states is
$$
\Omega(E=m \epsilon) = \frac{N!}{m! (N-m)!}
$$
Use the method of coincidences outlined in class to illustrate why computing the quantity $\Omega$ is so difficult. Consider a small number of particles $N = 50$, and then place $m$ of them in the excited state, and $N-m$ in the ground state. Coincidences here are meant that for any pair of particles, they are both in the same state. The ratio (total # pairs/# coincidences) is your estimate of $\Omega$.
Compute $\Omega$ for the low lying energy levels (small $m$). See how high in $m$ you can go before the result becomes unreliable.
Note: be careful about calculating factorial of large numbers.
!!!! Hint
One state of the system is a specific list $(1,0,1,0, \ldots ,0)$, which is $N$ elements long. In this state, the first particle is in the excited state, the second in the ground state, etc. Any state of fixed energy (fixed $m$) will have $m$ ones, and $N-m$ zeroes.
In order to estimate $\Omega$ you need to generate a large number of states (say 10,000 such lists, each with $N$ entries), and then determine how many coincidences occur. That is, a coincidence is the appearance of the same state (list) twice. Of course, the probability of such a coincidence is ridiculously small, $1/\Omega$.
!!!!Homeworks
* [[Homework 1|Homework 1 2020]]. Due on Tuesday, September 22, 2020.
* [[Homework 2|Homework 2 2020]]. Due on Wednesday, September 30, 2020.
* [[Homework 3|Homework 3 2020]]. Due on Wednesday, October 7, 2020.
* [[Homework 4|Homework 4 2020]]. Due on Wednesday, October 14, 2020.
* [[Homework 5|Homework 5 2020]]. Due on Wednesday, October 21, 2020.
!!!!Laboratory sessions
| ''8:00 a.m. - 8:50 a.m.'' | | |
| | Session A1 | https://umn.zoom.us/j/99254944724 |
| | Session B1 | https://umn.zoom.us/j/7695447358 |
| | Session C1 | https://umn.zoom.us/j/3810472354 |
| ''9:00 a.m. - 9:50 a.m'' |||
| | Session A2 | https://umn.zoom.us/j/91679894787 |
| | Session B2 | https://umn.zoom.us/j/7695447358 |
| | Session C2 | https://umn.zoom.us/j/3810472354 |
* [[Laboratory Session 1|LabSession 1 2020]]. September 10, 2020.
* [[Laboratory Session 2|LabSession 2 2020]]. September 17, 2020.
* [[Laboratory Session 3|LabSession 3 2020]]. September 24, 2020.
* [[Laboratory Session 4|LabSession 4 2020]]. October 1, 2020.
* [[Laboratory Session 5|LabSession 5 2020]]. October 8, 2020.
* [[Laboratory Session 6|LabSession 6 2020]]. October 15, 2020.
* [[Laboratory Session 7|LabSession 7 2020]]. October 22, 2020.
!!!! Mid term and Final exams
* Both center around a semester long computational project. [[There is a list of suggested projects|Exam Projects]]. A certain amount of discussion with the instructor will be necessary to help you select an appropriate project. It is also possible for you to suggest your own computational project in which you are interested (please, not your thesis research if you are a graduate student).
* The report in lieu of the Mid Term exam is due on October 31, 2020.
* The final report is due on December 17, 2020.
/***
|Name|ImageSizePlugin|
|Source|http://www.TiddlyTools.com/#ImageSizePlugin|
|Version|1.2.3|
|Author|Eric Shulman|
|License|http://www.TiddlyTools.com/#LegalStatements|
|~CoreVersion|2.1|
|Type|plugin|
|Description|adds support for resizing images|
This plugin adds optional syntax to scale an image to a specified width and height and/or interactively resize the image with the mouse.
!!!!!Usage
<<<
The extended image syntax is:
{{{
[img(w+,h+)[...][...]]
}}}
where ''(w,h)'' indicates the desired width and height (in CSS units, e.g., px, em, cm, in, or %). Use ''auto'' (or a blank value) for either dimension to scale that dimension proportionally (i.e., maintain the aspect ratio). You can also calculate a CSS value 'on-the-fly' by using a //javascript expression// enclosed between """{{""" and """}}""". Appending a plus sign (+) to a dimension enables interactive resizing in that dimension (by dragging the mouse inside the image). Use ~SHIFT-click to show the full-sized (un-scaled) image. Use ~CTRL-click to restore the starting size (either scaled or full-sized).
<<<
!!!!!Examples
<<<
{{{
[img(100px+,75px+)[images/meow2.jpg]]
}}}
[img(100px+,75px+)[images/meow2.jpg]]
{{{
[<img(34%+,+)[images/meow.gif]]
[<img(21% ,+)[images/meow.gif]]
[<img(13%+, )[images/meow.gif]]
[<img( 8%+, )[images/meow.gif]]
[<img( 5% , )[images/meow.gif]]
[<img( 3% , )[images/meow.gif]]
[<img( 2% , )[images/meow.gif]]
[img( 1%+,+)[images/meow.gif]]
}}}
[<img(34%+,+)[images/meow.gif]]
[<img(21% ,+)[images/meow.gif]]
[<img(13%+, )[images/meow.gif]]
[<img( 8%+, )[images/meow.gif]]
[<img( 5% , )[images/meow.gif]]
[<img( 3% , )[images/meow.gif]]
[<img( 2% , )[images/meow.gif]]
[img( 1%+,+)[images/meow.gif]]
{{tagClear{
}}}
<<<
!!!!!Revisions
<<<
2011.09.03 [1.2.3] bypass addStretchHandlers() if no '+' suffix is used (i.e., not resizable)
2010.07.24 [1.2.2] moved tip/dragtip text to config.formatterHelpers.imageSize object to enable customization
2009.02.24 [1.2.1] cleanup width/height regexp, use '+' suffix for resizing
2009.02.22 [1.2.0] added stretchable images
2008.01.19 [1.1.0] added evaluated width/height values
2008.01.18 [1.0.1] regexp for "(width,height)" now passes all CSS values to browser for validation
2008.01.17 [1.0.0] initial release
<<<
!!!!!Code
***/
//{{{
version.extensions.ImageSizePlugin= {major: 1, minor: 2, revision: 3, date: new Date(2011,9,3)};
//}}}
//{{{
var f=config.formatters[config.formatters.findByField("name","image")];
f.match="\\[[<>]?[Ii][Mm][Gg](?:\\([^,]*,[^\\)]*\\))?\\[";
f.lookaheadRegExp=/\[([<]?)(>?)[Ii][Mm][Gg](?:\(([^,]*),([^\)]*)\))?\[(?:([^\|\]]+)\|)?([^\[\]\|]+)\](?:\[([^\]]*)\])?\]/mg;
f.handler=function(w) {
this.lookaheadRegExp.lastIndex = w.matchStart;
var lookaheadMatch = this.lookaheadRegExp.exec(w.source)
if(lookaheadMatch && lookaheadMatch.index == w.matchStart) {
var floatLeft=lookaheadMatch[1];
var floatRight=lookaheadMatch[2];
var width=lookaheadMatch[3];
var height=lookaheadMatch[4];
var tooltip=lookaheadMatch[5];
var src=lookaheadMatch[6];
var link=lookaheadMatch[7];
// Simple bracketted link
var e = w.output;
if(link) { // LINKED IMAGE
if (config.formatterHelpers.isExternalLink(link)) {
if (config.macros.attach && config.macros.attach.isAttachment(link)) {
// see [[AttachFilePluginFormatters]]
e = createExternalLink(w.output,link);
e.href=config.macros.attach.getAttachment(link);
e.title = config.macros.attach.linkTooltip + link;
} else
e = createExternalLink(w.output,link);
} else
e = createTiddlyLink(w.output,link,false,null,w.isStatic);
addClass(e,"imageLink");
}
var img = createTiddlyElement(e,"img");
if(floatLeft) img.align="left"; else if(floatRight) img.align="right";
if(width||height) {
var x=width.trim(); var y=height.trim();
var stretchW=(x.substr(x.length-1,1)=='+'); if (stretchW) x=x.substr(0,x.length-1);
var stretchH=(y.substr(y.length-1,1)=='+'); if (stretchH) y=y.substr(0,y.length-1);
if (x.substr(0,2)=="{{")
{ try{x=eval(x.substr(2,x.length-4))} catch(e){displayMessage(e.description||e.toString())} }
if (y.substr(0,2)=="{{")
{ try{y=eval(y.substr(2,y.length-4))} catch(e){displayMessage(e.description||e.toString())} }
img.style.width=x.trim(); img.style.height=y.trim();
if (stretchW||stretchH) config.formatterHelpers.addStretchHandlers(img,stretchW,stretchH);
}
if(tooltip) img.title = tooltip;
// GET IMAGE SOURCE
if (config.macros.attach && config.macros.attach.isAttachment(src))
src=config.macros.attach.getAttachment(src); // see [[AttachFilePluginFormatters]]
else if (config.formatterHelpers.resolvePath) { // see [[ImagePathPlugin]]
if (config.browser.isIE || config.browser.isSafari) {
img.onerror=(function(){
this.src=config.formatterHelpers.resolvePath(this.src,false);
return false;
});
} else
src=config.formatterHelpers.resolvePath(src,true);
}
img.src=src;
w.nextMatch = this.lookaheadRegExp.lastIndex;
}
}
config.formatterHelpers.imageSize={
tip: 'SHIFT-CLICK=show full size, CTRL-CLICK=restore initial size',
dragtip: 'DRAG=stretch/shrink, '
}
config.formatterHelpers.addStretchHandlers=function(e,stretchW,stretchH) {
e.title=((stretchW||stretchH)?this.imageSize.dragtip:'')+this.imageSize.tip;
e.statusMsg='width=%0, height=%1';
e.style.cursor='move';
e.originalW=e.style.width;
e.originalH=e.style.height;
e.minW=Math.max(e.offsetWidth/20,10);
e.minH=Math.max(e.offsetHeight/20,10);
e.stretchW=stretchW;
e.stretchH=stretchH;
e.onmousedown=function(ev) { var ev=ev||window.event;
this.sizing=true;
this.startX=!config.browser.isIE?ev.pageX:(ev.clientX+findScrollX());
this.startY=!config.browser.isIE?ev.pageY:(ev.clientY+findScrollY());
this.startW=this.offsetWidth;
this.startH=this.offsetHeight;
return false;
};
e.onmousemove=function(ev) { var ev=ev||window.event;
if (this.sizing) {
var s=this.style;
var currX=!config.browser.isIE?ev.pageX:(ev.clientX+findScrollX());
var currY=!config.browser.isIE?ev.pageY:(ev.clientY+findScrollY());
var newW=(currX-this.offsetLeft)/(this.startX-this.offsetLeft)*this.startW;
var newH=(currY-this.offsetTop )/(this.startY-this.offsetTop )*this.startH;
if (this.stretchW) s.width =Math.floor(Math.max(newW,this.minW))+'px';
if (this.stretchH) s.height=Math.floor(Math.max(newH,this.minH))+'px';
clearMessage(); displayMessage(this.statusMsg.format([s.width,s.height]));
}
return false;
};
e.onmouseup=function(ev) { var ev=ev||window.event;
if (ev.shiftKey) { this.style.width=this.style.height=''; }
if (ev.ctrlKey) { this.style.width=this.originalW; this.style.height=this.originalH; }
this.sizing=false;
clearMessage();
return false;
};
e.onmouseout=function(ev) { var ev=ev||window.event;
this.sizing=false;
clearMessage();
return false;
};
}
//}}}
[img(100%+,auto)[http://homepages.spa.umn.edu/~vinals/tspot_files/Kmeans1d.jpg]]
[img(100%+,auto)[http://homepages.spa.umn.edu/~vinals/tspot_files/Kmeans2d.jpg]]
# Try your login credentials, or create them if you have not done so.
# Write a simple code, compile it an run it that will estimate the largest simple precision (4 byte) and double precision (8 byte) floating point and integer numbers that the computer you are using can represent. Hint: Use a loop to compute 1*2*2*2 .... until the returned value is out of range.
# The Bessel function of order 1 is given by $$ J_{1}(x) = \frac{\sin x}{x^{2}} - \frac{\cos x}{x}$$ Given the Taylor series for sin and cos, $$ \sin (x) = x - \frac{x^{3}}{3!} + \frac{x^{5}}{5!} - ... \quad \cos (x) = 1 - \frac{x^{2}}{2!} + \frac{x^{4}}{4!} - ...$$Write a code that computes this function in $x \in [0,10]$ and plot the result.
# Try your x500 login credentials. Ideally you should have tried them before you come to the lab.
# Log in github.umn.edu with your x500 credentials. Inform the instructor by e-mail when you have done so.
# Write a simple code, compile it an run it that will estimate the largest simple precision (4 byte) and double precision (8 byte) floating point and integer numbers that the computer you are using can represent. Hint: Use a loop to compute 1*2*2*2 .... until the returned value is out of range.
# The Bessel function of order 1 is given by $$ J_{1}(x) = \frac{\sin x}{x^{2}} - \frac{\cos x}{x}$$ Given the Taylor series for sin and cos, $$ \sin (x) = x - \frac{x^{3}}{3!} + \frac{x^{5}}{5!} - ... \quad \cos (x) = 1 - \frac{x^{2}}{2!} + \frac{x^{4}}{4!} - ...$$Write a code that uses the Taylor series to compute the Bessel function in $x \in [0,10]$ to a reasonable approximation and plot the result.
{{justifyright{
''Thursday, September 6, 2018''
}}}
# Try your x500 login credentials. Ideally you should have tried them before you come to the lab.
# Log in github.umn.edu with your x500 credentials. Inform the instructor by e-mail when you have done so.
# Write a simple code, compile it an run it that will estimate the largest simple precision (4 byte) and double precision (8 byte) floating point and integer numbers that the computer you are using can represent. Hint: Use a loop to compute 1*2*2*2 .... until the returned value is out of range.
# The Bessel function of order 1 is given by $$ J_{1}(x) = \frac{\sin x}{x^{2}} - \frac{\cos x}{x}$$ Given the Taylor series for sin and cos, $$ \sin (x) = x - \frac{x^{3}}{3!} + \frac{x^{5}}{5!} - ... \quad \cos (x) = 1 - \frac{x^{2}}{2!} + \frac{x^{4}}{4!} - ...$$Write a code that uses the Taylor series to compute the Bessel function in $x \in [0,10]$ to a reasonable approximation and plot the result.
Compilation instructions:
* ''FORTRAN'' (type man gfortran for information and options. Information [[about Fortran file extensions|http://fortranwiki.org/fortran/show/File+extensions]]).
** gfortran filename.f90
* ''C/C++'' (type man g++ or man gcc for information and options)
** g++ filename.c or gcc filename.c -lm
{{justifyright{
''Thursday, September 10, 2020''
}}}
# Try your UM ID login credentials on the ~CSELabs systems. Remotely, the easiest way is to access https://vole.cse.umn.edu. Ideally you have done this prior to this session. Note that this will open a Unix session.
# Log in github.umn.edu with your UM ID. Log out. Inform the instructor by e-mail when you have done so.
# Two small programming exercises to get you acquainted with file editors, compilers, and execution environments. You can create and run the code on the ~CSELabs systems, or on your own personal computer. In this case, you will need to install the appropriate compilers.
## Write a simple code, compile it an run it that will estimate the largest simple precision (4 byte) and double precision (8 byte) floating point and integer numbers that the computer you are using can represent. Hint: Use a loop to compute 1*2*2*2 .... until the returned value is out of range.
## The Bessel function of order 1 is given by $$ J_{1}(x) = \frac{\sin x}{x^{2}} - \frac{\cos x}{x}$$ Given the Taylor series for sin and cos, $$ \sin (x) = x - \frac{x^{3}}{3!} + \frac{x^{5}}{5!} - ... \quad \cos (x) = 1 - \frac{x^{2}}{2!} + \frac{x^{4}}{4!} - ...$$Write a code that uses the Taylor series to compute the Bessel function in $x \in [0,10]$ to a reasonable approximation and plot the result.
Compilation instructions:
* ''FORTRAN'' (type man gfortran for information and options. Information [[about Fortran file extensions|http://fortranwiki.org/fortran/show/File+extensions]]).
** gfortran filename.f90
* ''C/C++'' (type man g++ or man gcc for information and options)
** g++ filename.c or gcc filename.c -lm
!!Random numbers and their distributions
So called pseudo-random numbers can be generated on a computer. There are multiple routines that can accomplih this to varying degrees of sophistication. We will use a simple intrinsic function on Linux computers. The function call (both FORTRAN and C)
{{{
x=RAND()
}}}
assigns to x a random number uniformly distributed between $[0,1)$ (any number in that interval with equal probability).
Because the numbers are not truly random (rather they are a predictable sequence), a sequence has to be initialized with an integer number called 'seed'. For example
{{{
INTEGER SEED
SEED = 123456
CALL SRAND(SEED)
...
x=RAND()
}}}
Every time the same seed is used, the same sequence of random numbers results. A sequence of random numbers is obtained simply by repeating the function call RAND().
!!!!1. Uniform random number statistics
Just to get started, generate 100,000 random numbers and compute the average and standard deviation of the sample. You should be able to predict the values analytically and compare with the results.
!!!!2. Histogram
We will create a histogram of the sequence of random numbers and compare it to the expected probability (uniform). In order to create a histogram, you can use the following procedure
{{{
INTEGER :: hist(100)
REAL(8) :: rhist(100)
DO i=1,100
hist(i)=0
END DO
dx=(1.0-0.0)/100
CALL SRAND(3452345)
DO i=1,100000
index=RAND()/dx
hist(index)=hist(index)+1
END DO
DO i=1,100
rhist(i)=hist(i)/(dx*100000)
WRITE(27,*) i*dx,rhist(i)
END DO
}}}
Plot the result and compare the histogram to the expected probability distribution of the random numbers.
!!!!3. Gaussianly distributed random numbers
Gaussianly distributed random numbers are generated in pairs (we will see this in class). Given two random numbers $r_{1}$ and $r_{2}$ uniformly distributed in [0,1), the following two numbers
$$x_{1} = \sqrt{-2 \ln r_{1}} \; \cos (2 \pi r_{2}) \quad\quad x_{2} = \sqrt{-2 \ln r_{1}} \; \sin (2 \pi r_{2})$$
are two independent random numbers, Gaussianly distributed ov zero average and variance equal to 1.
Repeat the process of (2) to generate the histogram of the resulting sequence.
!!!!4. Correlations
The sequences above involve independent random numbers and therefore there are no correlations. We will examine this first. Assume that you have generated the sequence of random numbers $\{ x_{i} \}$. The correlation between, say, $x_{i}$ and $x_{i+1}$ is defined
$$C_{i\;i+1} = \langle x_{i} x_{i+1} \rangle - \langle x_{i} \rangle \langle x_{i+1} \rangle$$
where average over $N$ points are defined as
$$\langle ... \rangle \approx \frac{1}{N} \sum_{i=1}^{N} (...)_{i}$$
Calculate $C_{i \; i+1}$ for the Gaussian sequence above.
Correlation functions are a generalization in that one defines
$$C(j) = \langle x_{i} x_{i+j} \rangle - \langle x_{i} \rangle \langle x_{i+j} \rangle$$
Note that the correlation only depends on $j$, the "distance" between the two points. This is normally the case in Physics, but it is not always the case.
!!!!5. More correlations
Consider the following sequence
$$x_{i+1} = x_{i} - x_{i-1} + r_{i}$$
where $r_{i}$ is a random number uniformly distributed in [0,1). Specify some initial condition $x_{0}$ and $x_{1}$. Run a long sequence and calculate the correlation function $C(j)$ for $j$ between 1 and 10.
* [[Presentation|http://homepages.spa.umn.edu/~vinals/tspot_files/Lab10.pdf]] for this session posted.
* We will next cover some [[Programming Tips]].
* The final code has been added to github.umn.edu.
{{justifyright{
''Thursday, November 8, 2018''
}}}
!!Random numbers and their distributions
So called pseudo-random numbers can be generated on a computer. There are multiple routines that can accomplish this to varying degrees of sophistication. We will use a simple intrinsic function on Linux computers. A function call (both in FORTRAN and C++) assigns a random number uniformly distributed between $[0,1)$ (any number in that interval with equal probability). Details on the syntax can be found [[here|http://homepages.spa.umn.edu/~vinals/classes/phys4041/phys4041_notes.pdf]].
Go over the sequence of steps indicated below as time allows during the laboratory session.
!!!!1. Uniform random number statistics
Just to get started, generate 100,000 random numbers and compute the average and standard deviation of the sample. You should be able to predict the values analytically and compare with the results.
!!!!2. Histogram
We will create a histogram of the sequence of random numbers and compare it to the expected probability (uniform). In order to create a histogram, you can use the following procedure (e.g., in FORTRAN. Similarly in C++)
{{{
INTEGER :: hist(100)
REAL(8) :: rhist(100)
DO i=1,100
hist(i)=0
END DO
dx=(1.0-0.0)/100
CALL SRAND(3452345)
DO i=1,100000
index=RAND()/dx + 1
hist(index)=hist(index)+1
END DO
DO i=1,100
rhist(i)=hist(i)/(dx*100000)
WRITE(27,*) i*dx,rhist(i)
END DO
}}}
Plot the result and compare the histogram to the expected probability distribution of the random numbers.
!!!!3. Gaussianly distributed random numbers
Gaussianly distributed random numbers are generated in pairs (we will see this in class). Given two random numbers $r_{1}$ and $r_{2}$ uniformly distributed in [0,1), the following two numbers
$$x_{1} = \sqrt{-2 \ln r_{1}} \; \cos (2 \pi r_{2}) \quad\quad x_{2} = \sqrt{-2 \ln r_{1}} \; \sin (2 \pi r_{2})$$
are two independent random numbers, Gaussianly distributed of zero average and variance equal to 1.
Repeat the process of (2) to generate the histogram of the resulting sequence.
!!!!4. Correlations
The sequences above involve independent random numbers and therefore there are no correlations. We will examine this first. Assume that you have generated the sequence of random numbers $\{ x_{i} \}$. The correlation between, say, $x_{i}$ and $x_{i+1}$ is defined
$$C_{i\;i+1} = \langle x_{i} x_{i+1} \rangle - \langle x_{i} \rangle \langle x_{i+1} \rangle$$
where average over $N$ points are defined as
$$\langle ... \rangle \approx \frac{1}{N} \sum_{i=1}^{N} (...)_{i}$$
Calculate $C_{i \; i+1}$ for the Gaussian sequence above.
Correlation functions are a generalization in that one defines
$$C(j) = \langle x_{i} x_{i+j} \rangle - \langle x_{i} \rangle \langle x_{i+j} \rangle$$
Note that the correlation only depends on $j$, the "distance" between the two points. This is normally the case in Physics, but it is not always the case.
!!!!5. More correlations
Consider the following sequence
$$x_{i+1} = x_{i} - x_{i-1} + r_{i}$$
where $r_{i}$ is a random number uniformly distributed in [0,1). Specify some initial condition $x_{0}$ and $x_{1}$. Run a long sequence and calculate the correlation function $C(j)$ for $j$ between 1 and 10.
!!!Correlation Functions
!!!!!Independent and Gaussianly distributed random numbers
Gaussianly distributed random numbers are generated in pairs (as seen in class). Given two random numbers $r_{1}$ and $r_{2}$ uniformly distributed in [0,1), the following two numbers
$$x_{1} = \sqrt{-2 \ln r_{1}} \; \cos (2 \pi r_{2}) \quad\quad x_{2} = \sqrt{-2 \ln r_{1}} \; \sin (2 \pi r_{2})$$
are two independent random numbers, Gaussianly distributed ov zero average and variance equal to 1. We wish to check both the Gaussianity and the independence.
In order to check for Gaussianity, compute $\langle x \rangle$, $\langle x^{2} \rangle$,$\langle x^{3} \rangle$ and $\langle x^{4} \rangle$ where the averages are defined over the sequence:
$$\langle ... \rangle \approx \frac{1}{N} \sum_{i=1}^{N} (...)_{i}$$
Show that odd moments are zero, that $\langle x^{2} \rangle = \sigma^{2} = 1$ (the variance of the sequence), and that $\langle x^{4} \rangle = 3 \langle x^{2} \rangle$ as appropriate for a Gaussian distribution.
In order to check for independence in a sequence of random numbers $\{ x_{i} \}$, we compute the correlation between, say, $x_{i}$ and $x_{i+1}$ defined by
$$C(1) = \langle x_{i} x_{i+1} \rangle - \langle x_{i} \rangle \langle x_{i+1} \rangle$$
Or more generally, we can define a correlation function
$$C(j) = \langle x_{i} x_{i+j} \rangle - \langle x_{i} \rangle \langle x_{i+j} \rangle$$
Calculate $C(j)$, $j=1, ..., 10$ for the sequence of Gaussian numbers above.
!!!!!Correlations in an oscillatory process
In general, stochastic processes generate states which are correlated. As an example, consider the sequence,
$$x_{i+1} = x_{i} - x_{i-1} + r_{i}$$
where $r_{i}$ is a sequence of independent random numbers uniformly distributed in [0,1). Specify some initial condition $x_{0}$ and $x_{1}$. This iteration produces a sequence of random numbers $\{ x_{i} \}$.
Run a long sequence and plot it. Next, calculate the correlation function $C(j)$ for $j$ between 1 and 10.
!!Random numbers and their distributions
So called pseudo-random numbers can be generated on a computer. There are multiple routines that can accomplih this to varying degrees of sophistication. We will use a simple intrinsic function on Linux computers. The function call (both FORTRAN and C)
{{{
x=RAND()
}}}
assigns to x a random number uniformly distributed between $[0,1)$ (any number in that interval with equal probability).
Because the numbers are not truly random (rather they are a predictable sequence), a sequence has to be initialized with an integer number called 'seed'. For example
{{{
INTEGER SEED
SEED = 123456
CALL SRAND(SEED)
...
x=RAND()
}}}
Every time the same seed is used, the same sequence of random numbers results. A sequence of random numbers is obtained simply by repeating the function call RAND().
Go over the sequence steps below as time allows during the laboratory session.
!!!!1. Uniform random number statistics
Just to get started, generate 100,000 random numbers and compute the average and standard deviation of the sample. You should be able to predict the values analytically and compare with the results.
!!!!2. Histogram
We will create a histogram of the sequence of random numbers and compare it to the expected probability (uniform). In order to create a histogram, you can use the following procedure
{{{
INTEGER :: hist(100)
REAL(8) :: rhist(100)
DO i=1,100
hist(i)=0
END DO
dx=(1.0-0.0)/100
CALL SRAND(3452345)
DO i=1,100000
index=RAND()/dx
hist(index)=hist(index)+1
END DO
DO i=1,100
rhist(i)=hist(i)/(dx*100000)
WRITE(27,*) i*dx,rhist(i)
END DO
}}}
Plot the result and compare the histogram to the expected probability distribution of the random numbers.
!!!!3. Gaussianly distributed random numbers
Gaussianly distributed random numbers are generated in pairs (we will see this in class). Given two random numbers $r_{1}$ and $r_{2}$ uniformly distributed in [0,1), the following two numbers
$$x_{1} = \sqrt{-2 \ln r_{1}} \; \cos (2 \pi r_{2}) \quad\quad x_{2} = \sqrt{-2 \ln r_{1}} \; \sin (2 \pi r_{2})$$
are two independent random numbers, Gaussianly distributed of zero average and variance equal to 1.
Repeat the process of (2) to generate the histogram of the resulting sequence.
!!!!4. Correlations
The sequences above involve independent random numbers and therefore there are no correlations. We will examine this first. Assume that you have generated the sequence of random numbers $\{ x_{i} \}$. The correlation between, say, $x_{i}$ and $x_{i+1}$ is defined
$$C_{i\;i+1} = \langle x_{i} x_{i+1} \rangle - \langle x_{i} \rangle \langle x_{i+1} \rangle$$
where average over $N$ points are defined as
$$\langle ... \rangle \approx \frac{1}{N} \sum_{i=1}^{N} (...)_{i}$$
Calculate $C_{i \; i+1}$ for the Gaussian sequence above.
Correlation functions are a generalization in that one defines
$$C(j) = \langle x_{i} x_{i+j} \rangle - \langle x_{i} \rangle \langle x_{i+j} \rangle$$
Note that the correlation only depends on $j$, the "distance" between the two points. This is normally the case in Physics, but it is not always the case.
!!!!5. More correlations
Consider the following sequence
$$x_{i+1} = x_{i} - x_{i-1} + r_{i}$$
where $r_{i}$ is a random number uniformly distributed in [0,1). Specify some initial condition $x_{0}$ and $x_{1}$. Run a long sequence and calculate the correlation function $C(j)$ for $j$ between 1 and 10.
{{justifyright{
''Thursday, November 15, 2018''
}}}
!!! Molecular Dynamics code for homework 8
Start the development of the code for homework 8: a Molecular Dynamics simulation of a two dimensional ~Lennard-Jones fluid. Depending on progress, go down the following list:
# Write a module that will compute all inter particle separations for an arbitrary configuration.
# Write another module that, given the interparticle separations, will compute the instantaneous force on each particle.
# Write a module that will implement the boundary conditions given.
# Write the main code that will conduct the time integration by using the Verlet algorithm.
# Write the analysis module in which, given a configuration, you will compute and output the desired quantities (energies, etc).
Continue debugging your Molecular Dynamics code.
Once it runs, experiment with:
* A small number of particles $N = 10$ in a small box $L=10$. Plot trajectories and examine the effect of boundary conditions. Hint: an easy way to plot trajectories is to write to a file with $(x,y)$ for all particles and a sequence of times and simply use gnuplot.
* Compare the trajectories between the boundary conditions given and a different condition in which particles are just reflected off boundaries.
* Examine the effect of temperature on the trajectories, as well as the values of the energy and pressure.
!!! Correlation function of the two dimensional ~Lennard-Jones system
Use the Molecular Dynamics code that you have written for the homework to compute the pairwise correlation function of the ~Lennard-Jones system for a range of temperatures. The correlation function is defined as a histogram of inter particle separations $$ g(\vec{r}) = \frac{V}{N} \langle \sum_{i,j} \delta(\vec{r} - \vec{r}_{ij}) \rangle
$$
where $\vec{r}_{ij}$ is the vector joining the positions of particles $i$ and $j$, and $j$ is different from $i$. The average is understood over time during the simulation. For simplicity, just focus on separations $r = \| \vec{r} \|$ and only compute the correlation function of separations $g(r)$.
Calculate and plot the correlation function $g(r)$ in the two following states: Fix the dimensionless temperature to $T = 0.45$ and examine two dimensionless densities $\rho = N/V = 0.855$ and $\rho = 0.786$.
{{justifyright{
''Thursday, November 29, 2018''
}}}
!!! Correlation function of the two dimensional ~Lennard-Jones system
Use the Molecular Dynamics code that you have written for the homework to compute the pairwise correlation function of the ~Lennard-Jones system for a range of temperatures. The correlation function is defined as a histogram of inter particle separations $$ g(\vec{r}) = \frac{S}{N^2} \langle \sum_{i, j} \delta(\vec{r} - \vec{r}_{ij}) \rangle
$$
where $\vec{r}_{ij}$ is the vector joining the positions of particles $i$ and $j$, ($j$ different from $i$), and $S$ is the total surface of your two dimensional system. The average is understood over time during the simulation. For simplicity, just focus on separations $r = \| \vec{r} \|$ and only compute the correlation function of separations $g(r)$.
Calculate and plot the correlation function $g(r)$ at the dimensionless temperature $T = 0.45$ and surface density $N/S = 0.786$.
!!Gillespie algorithm for a simple production-degradation reaction system
Consider the following two reactions giving production and degradation of X, $$\emptyset \;\; \overset{A}\rightarrow \;\; X \quad\quad X \;\; \overset{B}\rightarrow \;\; \emptyset$$
where production occurs at a constant rate $A$ and degradation at a constant rate $B$.
If $n$ is the number of molecules of $X$, the Master Equation that corresponds to these two reactions is,$$ \frac{\partial p_{n}(t)}{\partial t} = A(E^{-1}-1) p_{n}(t) + B (E-1)(np_{n}(t))$$
Identifying the propensitties of the two reactions as $a_{1} = A$ and $a_{2} = Bn$ write a code that implements the Gillespie algorithm to evolve some arbitrary initial condition to the steady state distribution.
In this case, the steady state distribution is known to be $$ p_{n} = \sqrt{ \frac{B}{2 \pi A}} e^{-\frac{B}{2A} \left( n - \frac{A}{B} \right)^{2}}$$
{{justifyright{
''Thursday, December 6, 2018''
}}}
!!! Stationary distribution of the ~Landau-Gilbert equation
The ~Landau-Gilbert equation governs the evolution of the three dimensional magnetization vector in an external field. Consider for simplicity a one dimensional chain $\vec{M}(x,t) = \left( M_{x}(x,t), M_{y}(x,t), M_{z}(x,t) \right) $. The equation is
\begin{equation}
\frac{1}{\gamma_{0}} \frac{\partial \vec{M}}{\partial t} + \eta \vec{M}
\times \frac{\partial \vec{M}}{\partial t} = \vec{M} \times \left[
\vec{H}_{0} + \frac{A}{M_{s}^{2}} \frac{\partial^{2} \vec{M}}{\partial x^{2}}
+ K_{e} \left( \vec{M}\cdot \hat{k} \right) \hat{k} \right]
\end{equation}
where $\gamma_{0}, \eta, \vec{H}_{0}, A, K_{e}$ and $M_{s} = \| \vec{M} \|$ are constants, assumed known, and $\hat{k}$ is the unit vector in the $z$ direction. The magnitude of the magnetization $M_{s}$ is constant (the vector $\vec{M}$ only changes orientation). In order to enforce this fact, the dynamic equation is supplemented with the constraint,
\begin{equation}
\frac{\partial \vec{M}}{\partial t} \cdot \vec{M} = 0
\end{equation}
Over time, the evolution of the magnetization given by the ~Landau-Gilbert equation leads to the minimum of the following energy,
\begin{equation}
E( \{ \vec{M}_{j} \} ) = \sum_{j=1}^{N-1} \left[ \frac{A}{2M_{s}^{2}} \left(
\vec{M}_{j+1} - \vec{M}_{j} \right)^{2} -\frac{1}{2} \vec{H}_{0} \cdot
\vec{M}_{j} - \frac{K_{e}}{2} \left(\vec{M}_{j} \cdot \hat{k} \right)^{2}
+ \lambda (\vec{M}_{j}\cdot\vec{M}_{j} - M_{s}^{2})^{2}\right]
\end{equation}
where $\lambda$ is a constant, and we have discretized space $x_{j} = j \Delta
x$, $j = 1, \ldots, N-1$.
Modify the Monte Carlo code that we have written in class to sample equilibrium configurations of $\{\vec{M}_{j} \}$ obeying the Boltzmann distribution associated with the energy $E( \{ \vec{M}_{j} \} )$ at temperature $T$ by using the Metropolis algorithm.
!Version control and respositories
We will use git as the version control, and the ~UofM github respository. Help on git can be found at https://www.atlassian.com/git
!!Create the respository where class code examples will be stored
* Point your browser to github.umn.edu and login using your X500 credentials.
* Clone the respository on your local machine (you can do this on as many machines as you wish. They will all syncronize). Go the directory in which you would like to have the respository created and type:
{{{
git clone https://github.umn.edu/jvinals/PHYS4041-2016.git
}}}
The subdirectory ~PHYS4041-2016 will be created which will be your local copy of the repository.
Although we will use Linux to conduct class activities, there are also desktop applications available for Mac and windows. You can also synchronize those operating systems if desired. The desktops can be found at https://desktop.github.com/
You can also create your own repository for your private work. To do so,
* log in github.umn.edu
* Click on green button "New repository".
* Give it a name (e.g., myrepo), and (typically) choose private.
* Select "Initialize with a README".
* On your local computer, go to the directory where you want the repository to reside and type
{{{
git clone https://github.umn.edu/[your x500 username]/myrepo.git
}}}
A new subdirectory myrepo will be created that can be syncronized with github.umn.edu.
!!Version control
!!!!!0. Synchronize
Synchronization between your local computer and the repository is simple. Once inside the directory that needs to be synchronized, you would type:
{{{
git pull
}}}
to synchronize the local files to those in the repository ("pull" from the repository).
!!!!!1. Add files that will be versioned to the repository
Use an editor to create the file test.f90. In order to turn on version control, you need to "add" the file to git. Type
{{{
git add test.f90
}}}
git maintains its own list of files subject to revision control. Other files can be created in the same directory that git will simply ignore.
In order to check the status of your files and versions, you can type
{{{
git status
}}}
In the case at hand, it will list a "new file", test.f90
!!!!!2. Committing changes
You can edit and manipulate your file at will. Once you are satisfied with the changes, you can "commit" the latest version to the local repository. git allows you to compare and restore all committed versions in the past. In order to commit, you would type
{{{
git commit
}}}
An editor window will open where you would type any informational message you want kept with the version committed. For example, "do loop fixed and print option added". Any message that will help you later understand what this version is. To finish, hit the Esc key followed by :wq (these are vi commands). Every file in the "add" list that has been changed has now been committed. You can type git status to see where you are.
!!!!!3. Branches
!!!!1. Use of make to maintain scientific code.
Write a makefile that will be used to compile the code you wrote for the homework to check the accuracy of the machine you are using.
Scientific code is usually comprised of a relatively large number of modules and links to external libraries. The command "make" in Linux simplifies code maintenance a great deal. I will list a couple of examples, here. A few more can be found [[here|http://linux.101hacks.com/unix/make/]].
The command make simply executes a file called (by default, this can be changed) makefile. The command simply interprets the statements contained in makefile. The simplest example concerns compiling the code that you produced to estimate the accuracy of the computer you are using. The task of compiling the code would be undertaken by make. Assuming your code was named accuracy.f, create a file makefile with:
{{{
accuracy: accuracy.f
[Tab] gfortran -o accuracy accuracy.f
}}}
Typing make at the command prompt, will compile the code accuracy.f ''only'' when modified.
The general format is:
{{{
Target : Dependencies
[Tab] command to execute if files in dependencies modified
}}}
Target is the final product of the make command. Target depends on other files (e.g. its own source code, or other external codes). If any of the dependencies is satisfied, then the command line following it (after the tab) will execute.
A more complex example in which two codes are part of the same task:
{{{
combined-file: file1.o file2.o
[Tab] gcc file1.o file2.o -o combined-file
file1.o : file1.c
[Tab] gcc -c file1.c -o file1.o
file2.o : file2.c
[Tab] gcc -c file2.c -o file2.o
}}}
make can parse variables, and do more complicated tasks. Here is a more complex example:
{{{
.SUFFIXES = .o .f90 .f
OPTION = OPTIM
COMP_FORTRAN = gfortran
COMP_C = gcc
OBJECTS = file1.o file2.o
final_code: final_code.f $(OBJECTS)
[Tab] $(COMP_FORTRAN) $(OPTION) -o final_code final_code.f $(OBJECTS)
.f.o:;
[Tab] $(COMP_FORTRAN) -c $.f
}}}
!!!!2. Git and version control
Modern day coding and code maintenance makes heavy use of version control software and repositories. For this class, we will use git as version control, and github.umn.edu as the repository. For this session we will only use it to download files.
# Open a browser, navigate to github.umn.edu and log on with your x500 credentials.
# @@color(maroon):Contact the instructor to receive access to the repository (after you have logged on at least once)@@.
# Search for the respository jvinals/phys4041.
# Click on the file to download.
# On the top right, click on Raw.
# Save the file from your browser to your computer.
You may also setup a local directory to copy an existing git repository.
# Open a terminal window on your machine.
# In your home directory (or create/navigate to one less cluttered), enter: {{{ git init }}}
# Enter: {{{ git clone https://github.umn.edu/jvinals/phys4041.git }}}
# You will need to enter your x500 credentials.
# Once downloaded, you will have a local copy of the repository in a directory named phys4041
# To download future updates to the repository, navigate to this directory and enter: {{{ git pull }}}
More info on git can be found in this tutorial: https://www.atlassian.com/git
!!!!3. Linking to external code or libraries
This is a simple test of using external code to accomplish the same task as in the homework. Either:
* From github.umn.edu/jvinals/phys4041 download the subroutines/functions dlamch.f and lsame.f (FORTRAN), dlamch.c, blaswrap.h (C).
* Write the main code that uses this subroutines/functions to obtain information about the accuracy of the computer. The source code contains help about the various parameters that the functions accept.
Or:
* ''Only on Physics systems'': You can link to the LAPACK package directly that contains those subroutines without any need to dowload. The entire package LAPACK (numerical linear algebra) already compiled is in the archive /usr/lib64/liblapack.so. If your main code were called dlamch_drive.f90, you would compile by typing:
{{{
gfortran dlamch_drive.f90 -o dlamch_drive /usr/lib64/liblapack.so
}}}
which is often abbreviated into
{{{
gfortran dlamch_drive.f90 -o dlamch_drive -L /usr/lib64/ -llapack
}}}
Create a makefile to automate the code build.
!!!!1. Use of make to maintain scientific code
Write a makefile that will be used to compile the code you wrote during the first Lab session to check the largest and smallest floating point numbers of the machine you are using.
Scientific code is usually comprised of a relatively large number of modules and links to external libraries. The command "make" in Linux simplifies code maintenance a great deal. I will list a couple of examples, here. A few more can be found [[here|http://linux.101hacks.com/unix/make/]].
The command make simply executes a file called (by default, this can be changed) makefile. The command simply interprets the statements contained in makefile. The simplest example concerns compiling the code mentioned above. The task of compiling the code would be undertaken by make. Assuming your code was named accuracy.f, create a file makefile with (FORTRAN)
{{{
accuracy: accuracy.f
[Tab] gfortran -o accuracy accuracy.f
}}}
or C:
{{{
accuracy: accuracy.c
[Tab] g++ -o accuracy accuracy.c
}}}
Typing make at the command prompt, will compile the code accuracy.f or accuracy.c ''only'' when modified, and create the executable named accuracy.
The general format of a makefile entry is:
{{{
Target: Dependencies
[Tab] command to execute if files in dependencies modified
[Tab] another command to be executed
}}}
Target is the final product of the make command. Target depends on other files (e.g. its own source code, or other external codes). If any of the dependencies is satisfied, then the command line following it (after the tab) will execute.
A more complex example in which two codes are part of the same task:
{{{
combined-file: file1.o file2.o
[Tab] g++ file1.o file2.o -o combined-file
file1.o: file1.c
[Tab] g++ -c file1.c -o file1.o
file2.o: file2.c
[Tab] g++ -c file2.c -o file2.o
}}}
make can parse variables, and do more complicated tasks. Here is a more complex example:
{{{
.SUFFIXES = .o .f90 .f
OPTION = OPTIM
COMP_FORTRAN = gfortran
COMP_C = gcc
OBJECTS = file1.o file2.o
final_code: final_code.f $(OBJECTS)
[Tab] $(COMP_FORTRAN) $(OPTION) -o final_code final_code.f $(OBJECTS)
.f.o:;
[Tab] $(COMP_FORTRAN) -c $.f
}}}
!!!!2. Linking to external code or libraries
In this case, you will learn how to use pieces of code (libraries) already written by others. As an example,
we will use already existing external code to find the accuracy of the machine that you are using.
Either:
* From github.umn.edu/jvinals/phys4041-2018 download the subroutines/functions dlamch.f and lsame.f (FORTRAN), dlamch.c, blaswrap.h (C). {{{git pull}}} should do it if you are set up correctly with github.
* Write the main code that uses this subroutines/functions to obtain information about the accuracy of the computer. The source code of the routines contains help about the various parameters that the functions accept.
Or (the best option):
* Alternatively, you can link to the existing LAPACK package directly that contains those subroutines without any need to download. The entire LAPACK package (numerical linear algebra) already compiled is in the archive (CSELABS) /usr/lib/liblapack.so for FORTRAN and /usr/lib/liblapacke.so (C), and (Physics) /usr/lib64/liblapack.so (FORTRAN) and /usr/lib64/liblapacke.so (C). For example, if your main code were called dlamch_drive.f90, you would compile it on a Physics department system by typing:
{{{
gfortran dlamch_drive.f90 -o dlamch_drive /usr/lib64/liblapack.so
}}}
which is often abbreviated into
{{{
gfortran dlamch_drive.f90 -o dlamch_drive -L /usr/lib64/ -llapack
}}}
In C, you would type something similar. In the CSELABS systems, you would type
{{{
gcc dlamch_drive.c -L /usr/lib -llapack
}}}
Create a makefile to automate the code build.
!!!1. Use of make to maintain scientific code
Goal: To write a "makefile" that will be used to compile the code you wrote during the first Lab session to check the largest and smallest floating point numbers of the machine you are using. Scientific code is usually comprised of a relatively large number of modules and links to external libraries. The command "make" in Linux simplifies code maintenance a great deal. I will list a couple of examples here. A few more can be found [[here|http://linux.101hacks.com/unix/make/]].
The command make simply executes a file called (by default, this can be changed) makefile. The command simply interprets the statements contained in makefile. The simplest example concerns compiling the code mentioned above. Assuming your code was named accuracy.f90 or accuracy.c, use a text editor to create a file makefile with the following two lines (FORTRAN)
{{{
accuracy: accuracy.f90
[Tab] gfortran -o accuracy accuracy.f90
}}}
or C:
{{{
accuracy: accuracy.c
[Tab] g++ -o accuracy accuracy.c
}}}
Typing make at the command prompt, will compile the code accuracy.f or accuracy.c ''only'' when modified, and create the executable named accuracy.
The general format of a makefile entry is:
{{{
Target: Dependencies
[Tab] command to execute if files in dependencies modified
[Tab] another command to be executed
}}}
Target is the final product of the make command. Target depends on other files (e.g. its own source code, or other external codes). If any of the dependencies is satisfied, then the command line following it (after the tab) will execute.
A more complex example in which two codes are part of the same task:
{{{
combined-file: file1.o file2.o
[Tab] g++ file1.o file2.o -o combined-file
file1.o: file1.c
[Tab] g++ -c file1.c -o file1.o
file2.o: file2.c
[Tab] g++ -c file2.c -o file2.o
}}}
make can parse variables, and do more complicated tasks. Here is a more complex example:
{{{
.SUFFIXES = .o .f90 .f
OPTION = OPTIM
COMP_FORTRAN = gfortran
COMP_C = gcc
OBJECTS = file1.o file2.o
final_code: final_code.f90 $(OBJECTS)
[Tab] $(COMP_FORTRAN) $(OPTION) -o final_code final_code.f90 $(OBJECTS)
.f90.o;
[TAB] $(COMP_FORTRAN) -c $.f90
.f.o:;
[Tab] $(COMP_FORTRAN) -c $.f
}}}
!!!2. Linking to external code or libraries
In this case, you will learn how to use pieces of code ("libraries") already written by others. As an example, we will use already existing external code to find the accuracy of the machine that you are using.
The most complicated part of the task is determining how to write your code so that it properly interfaces with the required library. All libraries have documentation, and you want to pay special attention to the arguments that are required and returned. In the case at hand:
* In C. Open with a text editor the file dlamch.c Line 54 describes the purpose of this library. Line 59 and following describes the arguments needed to interface with the library. Line 21 gives the entry point to the library, which indicates how you should call it from your code.
* In FORTRAN. Open with a text editor the file dlamch.f. The arguments are described in line 20 and following. How to call it can be determined from line 1.
Either:
* From github.umn.edu/jvinals/~PHYS4041-2018 download the subroutines/functions dlamch.f and lsame.f (FORTRAN), dlamch.c, lsame.c blaswrap.h (C) that are part of the LAPACK package. {{{git pull}}} should do it if you are set up correctly with github.
* Write the main code that uses this subroutines/functions to obtain information about the accuracy of the computer.
* To compile your code, e.g., dlam_drive.c or dlam_drive.f90, and produce the executable dlamch_drive type
** FORTRAN
{{{
gfortran dlamch_drive.f90 -o dlamch_drive dlamch.f lsame.f
}}}
** C
{{{
gcc dlamch_drive.c -o dlamch_drive dlamch.c lsame.c /usr/lib/x86_64-linux-gnu/libf2c.so -u MAIN__
}}}
Or:
* Alternatively, you can directly link to the LAPACK set of libraries that is already installed in your computer. The compiler will find and use the subroutines/functions that you need without any need to download them one by one. The entire LAPACK package (numerical linear algebra) already compiled is in the archive (CSELABS) /usr/lib/liblapack.so for FORTRAN and /usr/lib/liblapacke.so (C), and (Physics) /usr/lib64/liblapack.so (FORTRAN) and /usr/lib64/liblapacke.so (C). For example, if your main code were called dlamch_drive.f90, you would compile it on a Physics department system, for example, by typing:
{{{
gfortran dlamch_drive.f90 -o dlamch_drive /usr/lib64/liblapack.so
}}}
which is often abbreviated into
{{{
gfortran dlamch_drive.f90 -o dlamch_drive -L /usr/lib64/ -llapack
}}}
In C, you would type something similar. In the CSELABS systems, for example, you would type
{{{
gcc dlamch_drive.c -o dlamch_drive -L /usr/lib -llapack
}}}
Create a makefile to automate the code build.
{{justifyright{
''Thursday, September 17, 2020''
}}}
!!! Task 1. Linking to external code or libraries
In a large number of computing tasks, you will want to reuse code written by others, or "libraries" either publicly available on the web, or commercial. We will use a first example from LAPACK (the linear algebra package in FORTRAN and C/C++) or numpy.linalg in Python. You will write code to solve a generic linear system of $n$ equations $A x = b$. In particular, use your code to solve
\begin{eqnarray}
3 x_{1} + 2 x_{2} & = & 2 \\
2 x_{1} + 6 x_{2} & = & -8
\end{eqnarray}
The most complicated part of the task is determining how to write your code so that it properly interfaces with the required library. All libraries have documentation, and you want to pay special attention to the arguments that are required and returned.
!!!! Information about LAPACK and numpy.linalg
* General user guide (for reference, do not start here): http://www.netlib.org/lapack/lug/
* General information for the C interface for Lapack: http://www.netlib.org/lapack/lapacke.html
* Information about [[numpy.linalg|https://numpy.org/doc/stable/reference/routines.linalg.html]].
* We will be using the routine {{{dgesv}}}. Documentation (FORTRAN) for the routine can be found [[here|http://www.netlib.org/lapack/double/dgesv.f]]. If you use C/C++, call the function {{{dgesv_(&n,&nrhs,a,&lda,ipiv,b,&ldb,&info)}}} where the length of array a is $lda*n$ (note also the underscore). Documentation [[here|http://www.physics.utah.edu/~detar/phys6720/handouts/lapack.html]]. The arguments and usage information are the same in both languages.
!!!! Coding instructions - FORTRAN users
* The LAPACK package on CSE systems is at: /lib/x86_64-linux-gnu/lapack/liblapack.so . There is an alias at /usr/lib as given below.
* You would compile your code with with the libraries by typing:
{{{
gfortran your_code.f90 -o your_code_executable -L /usr/lib/ -llapack
}}}
!!!! Coding instructions - C/C++ users
* The LAPACKE package on CSE systems is at: /lib/x86_64-linux-gnu/liblapacke.so . There is an alias at /usr/lib as given below.
* The libraries have been written in C. If your code is in C++, you have to alert the compiler that you will be linking to a library written in C. The necessary statement in your code is:
{{{
//declare function(s) from LAPACK
extern "C" void dgesv_(long*, long*, double*, long*, long*, double*, long*, long*);
}}}
* You would then compile your code with:
{{{
g++ your_code.c -o your_code_executable -L /usr/lib -llapack
}}}
!!!! Coding instructions - Python users
* The ~NumPy libraries are included in the Python installations (both 2.7 and 3.8) on the CSE systems. The typical convention is to
{{{
import numpy as np
}}}
at the start of the script; linalg functions would the be called as {{{np.linalg.f()}}}.
* ~NumPy's implementation runs on BLAS and LAPACK under the hood; the function {{{numpy.linalg.solve}}} is in fact implemented using {{{dgesv}}} and is thus roughly equivalent (though it provides somewhat less information due to lack of 'output parameters' such as {{{ipiv, info}}}).
!!! Task 2. Spherical Bessel function approximation (if you finish Task 1 early)
Spherical Bessel functions are trancendental functions that can be obtained by summing a series. However, there are recursion formulae that make the computation more efficient, although, recursion can introduce its own numerical problems. Use the following recursion formulae to write a code that computes the Bessel function $J_{l}$ with $l = 0, 1, 2, 3, 4$ and plot the results for $x \in (0,12)$. The recursion formulae ("upwards" or "downwards") are \begin{eqnarray}
J_{l+1}(x) & = & \frac{2l+1}{x} J_{l}(x) - J_{l-1}(x) \\
J_{l-1}(x) & = & \frac{2l+1}{x} J_{l}(x) - J_{l+1}(x) \\
J_{0}(x) = \frac{\sin x}{x} & \; & J_{1}(x) = \frac{\sin x - x \cos x}{x^{2}}
\end{eqnarray} Using the following table to help evaluate the accuracy of your approximations,
|! x |! $J_{3}$ |! $J_{5}$ |! $J_{8}$ |
| 0.10 | 9. 518519719 10$^{-6}$ | 9.616310231 10$^{-10}$ | 2.901200102 10$^{-16}$ |
| 1.00 | 9.006581118 10$^{-3}$ | 9.256115862 10$^{-5}$ | 2.826498802 10$^{-8}$ |
| 10.0 | -3.949584498 10$^{-2}$ | -5.553451162 10$^{-1}$ | 1.255780236 10$^{-1}$ |
** Compute the spherical Bessel functions required by using the upward recursion.
** Compare your solution to the table. Argue how errors might have propagated and what effect this may have on the accuracy of the solution.
** Chose an arbitrary starting function $J_{L}(x)$ with large $L$ and use the iteration downward. Compare the results to the upward iteration for the spherical Bessel functions asked above. Note that the downward recursion relation is unchanged if you multiply both sides by an arbitrary function $f(x)$. As a consequence, for an arbitrary starting function, the sequence obtained is defined up to this factor. Let $J'_{0}(x)$ be the zero-th order obtained, and $J_{0}(x)$ the exact function as given. Compute the ratio $J_{0}(x)/J'_{0}(x)$ and scale ''all'' functions $J_{l}$ obtained by this ratio.
!!!!Combination of spline finding and ~Newton-Raphson method
In problem 2 of Homework three, you wrote a code that uses spline interpolation to find, given the data, the volume $v$ at which the interpolation gives $p=3.25$.
Recast the problem as finding the solution to $p(v) - 3.25 = 0$. Use the ~Newton-Raphson algorithm to accomplish this ''with'' the spline interpolant as the function.
Make sure you write the appropriate makefile to compile this set of codes.
!!!!1. Bessel function approximation
Bessel functions are trancendental functions that can be obtained by summing a series. However, there are recursion formulae that make the computation more efficient, although, as seen in class, recursion can introduce its own numerical problems. Use the following recursion formulae to write a code that computes the Bessel function $J_{l}$ with $l = 0, 1, 2, 3, 4$ and plot the results for $x \in (0,12)$. The recursion formulae ("upwards" or "downwards") are \begin{eqnarray}
J_{l+1}(x) & = & \frac{2l+1}{x} J_{l}(x) - J_{l-1}(x) \\
J_{l-1}(x) & = & \frac{2l+1}{x} J_{l}(x) - J_{l+1}(x) \\
J_{0}(x) = \frac{\sin x}{x} & \; & J_{1}(x) = \frac{\sin x - x \cos x}{x^{2}}
\end{eqnarray} Using the following table to help evaluate the accuracy of your approximations,
|! x |! $J_{3}$ |! $J_{5}$ |! $J_{8}$ |
| 0.10 | 9. 518519719 10$^{-6}$ | 9.616310231 10$^{-10}$ | 2.901200102 10$^{-16}$ |
| 1.00 | 9.006581118 10$^{-3}$ | 9.256115862 10$^{-5}$ | 2.826498802 10$^{-8}$ |
| 10.0 | -3.949584498 10$^{-2}$ | -5.553451162 10$^{-1}$ | 1.255780236 10$^{-1}$ |
** Compute the Bessel functions required by using the upward recursion.
** Compare your solution to the table. Argue how errors might have propagated and what effect this may have on the accuracy of the solution.
** Chose an arbitrary starting function $J_{L}(x)$ with large $L$ and use the iteration downward. Compare the results to the upward iteration for the Bessel functions asked above. Note that the downward recursion relation is unchanged if you multiply both sides by an arbitrary function $f(x)$. As a consequence, for an arbitrary starting function, the sequence obtained is defined up to this factor. Let $J'_{0}(x)$ be the zero-th order obtained, and $J_{0}(x)$ the exact function as given. Compute the ratio $J_{0}(x)/J'_{0}(x)$ and scale ''all'' functions $J_{l}$ obtained by this ratio.
{{justifyright{
''Thursday, September 20, 2018''
}}}
A simple procedure to calculate numerically the //sin// function is through its Taylor series
$$
\sin(x) \approx \sin_{N}(x) = \sum_{l=0}^{l=N} (-1)^{l} \frac{x^{2l+1}}{(2l+1)!} = \sum_{l=0}^{N} A_{l}
$$
Write a code that that approximates $\sin(x)$ by $\sin_{13}(x)$ for 100,000 points uniformly distributed in $[ -\pi,\pi ]$, that is, by the first 14 terms in the Taylor series expansion.
As was the case in the homework, this is an alternating series, and it can be evaluated more efficiently through recursion. In particular, we have the result
$$
A_{l} = - \frac{x^{2}}{2l(2l+1)} A_{l-1} \quad {\rm and} \; A_{0} = x
$$
# Compute the same values of $sin(x)$ by this recursion, and compare them to what you obtained by summing the series directly (e.g., by plotting both).
# Compare the execution times of the two methods, and compare them to using directly the intrinsic function $sin(x)$ provided by your computer.
!!!!! Computing elapsed execution time in FORTRAN
{{{
CALL cpu_time(t0) ! Clock starts
! ... tasks to be timed ...
CALL cpu_time(t1)
WRITE (6,*) 'Elapsed time (in seconds): ',t1-t0
}}}
!!!!! Computing elapsed execution time in C++
{{{
#include <time.h>
clock_t time; // start clock time
double diff; // interval time in sec
.....
time = clock(); // start the clock
// Include statements that you want to time
diff = (double)(difftime(clock(), time))/CLOCKS_PER_SEC; // elapsed time since start
}}}
{{justifyright{
''Thursday, September 24, 2020''
}}}
!!! Task 1. (If you did not get to this last week): Spherical Bessel function approximation.
Spherical Bessel functions are trancendental functions that can be obtained by summing a series. However, there are recursion formulae that make the computation more efficient, although, recursion can introduce its own numerical problems. Use the following recursion formulae to write a code that computes the Bessel function $J_{l}$ with $l = 0, 1, 2, 3, 4$ and plot the results for $x \in (0,12)$. The recursion formulae ("upwards" or "downwards") are \begin{eqnarray}
J_{l+1}(x) & = & \frac{2l+1}{x} J_{l}(x) - J_{l-1}(x) \\
J_{l-1}(x) & = & \frac{2l+1}{x} J_{l}(x) - J_{l+1}(x) \\
J_{0}(x) = \frac{\sin x}{x} & \; & J_{1}(x) = \frac{\sin x - x \cos x}{x^{2}}
\end{eqnarray} Using the following table to help evaluate the accuracy of your approximations,
|! x |! $J_{3}$ |! $J_{5}$ |! $J_{8}$ |
| 0.10 | 9. 518519719 10$^{-6}$ | 9.616310231 10$^{-10}$ | 2.901200102 10$^{-16}$ |
| 1.00 | 9.006581118 10$^{-3}$ | 9.256115862 10$^{-5}$ | 2.826498802 10$^{-8}$ |
| 10.0 | -3.949584498 10$^{-2}$ | -5.553451162 10$^{-2}$ | 1.255780236 10$^{-1}$ |
** Compute the spherical Bessel functions required by using the upward recursion.
** Compare your solution to the table. Argue how errors might have propagated and what effect this may have on the accuracy of the solution.
** Chose two arbitrary starting functions $J_{L}(x)$ and $J_{L-1}(x)$ with large $L$ and use the iteration downward. Compare the results to the upward iteration for the spherical Bessel functions asked above. Note that the downward recursion relation is unchanged if you multiply both sides by an arbitrary function $f(x)$. As a consequence, for an arbitrary starting function, the sequence obtained is defined up to this factor. Let $J'_{0}(x)$ be the zero-th order obtained, and $J_{0}(x)$ the exact function as given. Compute the ratio $J_{0}(x)/J'_{0}(x)$ and scale ''all'' functions $J_{l}$ obtained by this ratio.
!!! Task 1 - new.
A simple procedure to calculate numerically the //sin// function is through its Taylor series
$$
\sin(x) \approx \sin_{N}(x) = \sum_{l=0}^{l=N} (-1)^{l} \frac{x^{2l+1}}{(2l+1)!} = \sum_{l=0}^{N} A_{l}
$$
Write a code that that approximates $\sin(x)$ by $\sin_{13}(x)$ for 100,000 points uniformly distributed in $[ -\pi,\pi ]$, that is, by the first 14 terms in the Taylor series expansion.
As was the case in the homework, this is an alternating series, and it can be evaluated more efficiently through recursion. In particular, we have the result
$$
A_{l} = - \frac{x^{2}}{2l(2l+1)} A_{l-1} \quad {\rm and} \; A_{0} = x
$$
# Compute the same values of $\sin(x)$ by this recursion, and compare them to what you obtained by summing the series directly (e.g., by plotting both).
# Compare the execution times of the two methods, and compare them to using directly the intrinsic function $\sin(x)$ provided by your computer.
!!!!! Computing elapsed execution time in FORTRAN
{{{
CALL cpu_time(t0) ! Clock starts
! ... tasks to be timed ...
CALL cpu_time(t1)
WRITE (6,*) 'Elapsed time (in seconds): ',t1-t0
}}}
!!!!! Computing elapsed execution time in C++
{{{
#include <time.h>
clock_t time; // start clock time
double diff; // interval time in sec
.....
time = clock(); // start the clock
// Include statements that you want to time
diff = (double)(difftime(clock(), time))/CLOCKS_PER_SEC; // elapsed time since start
}}}
!!!!! Computing elapased execution time in Python
{{{
from timeit import default_timer
start = default_timer();
# ... tasks to be timed ...
end = default_timer();
seconds_elapsed = end - start
}}}
You can also use the [[timeit module|https://docs.python.org/3/library/timeit.html#examples]] to time multiple iterations of the task in question in one go for better accuracy.
!!! Task 2 (if time allows). Physical (not numerical) instability
Imagine some population that starts out with $x_{0}$ individuals. The population grows with a growth rate $\mu$, but as the number becomes large, there is a limiting factor to growth (food, etc.). In discrete time steps, $n$, the governing equation for the evolution of the population is
$$
x_{n+1} = \mu x_{n}(1-x_{n})
$$
# Find the fixed points $x^{*}$ of this iteration ($\lim_{n\rightarrow \infty} x_{n} = x^{*}$).
# Iterate this equation for the following cases, and examine the temporal behavior: $\mu < 1$, $\mu = 2.9$, and $\mu = 3.8$.
!!Solution of a system of equations by using LAPACK
Lapack is an open source, widely used set of routines to accomplish a wide range of numerical linear algebra tasks. We will use it here to solve a linear system of the type $Ax =b$.
Information about Lapack:
* User Guide: http://www.netlib.org/lapack/lug/
* We will be using the routine dgesv. Detailed information at: http://www.netlib.org/lapack/explore-html/d8/d72/dgesv_8f.html#a5ee879032a8365897c3ba91e3dc8d512 This is the direct URL for this routine. A useful interface for specific routine usage is: http://www.netlib.org/lapack/explore-html/index.html
* C interface for Lapack: http://www.netlib.org/lapack/lapacke.html
!!!!Coding instructions - FORTRAN users
Computers in the Physics department (not CSELAB) have Lapack installed. You will have to ssh to a Physics host. Then when you compile, you need to link the Lapack libraries. All libraries can be linked at compilation time as
{{{
gfortran your_code.f90 -o your_code_executable -L /usr/lib64/ -llapack
}}}
!!!!Coding instructions - C users
Neither Physics computers nor CSELAB have the C version. We have uploaded the libraries into the class repository. The libraries are in liblapacke.a . Header files are also in the respository. Navigate to the directory in which you created the git repository for this class, and update it:
{{{
git pull
}}}
You would then link your code as
{{{
gcc/g++ your_code.c -o your_code_executable -llapack -lblas liblapacke.a -I include/
}}}
Note that the statement
{{{
#include "lapacke.h"
}}}
is needed in your code your_code.c. Note also that the function name is ~LAPACKE_dgesv, with the same calling arguments as the FORTRAN [[version|http://www.netlib.org/lapack/explore-html/d8/d72/dgesv_8f.html#a5ee879032a8365897c3ba91e3dc8d512]].
!!!!1. Solution of a simple system of linear equations
Write a code (and makefile) that using the Lapack routine dgesv solves the following system of equations:
\begin{eqnarray}
3 x_{1} + 2 x_{2} & = & 2 \\
2 x_{1} + 6 x_{2} & = & -8
\end{eqnarray}
!!!!2. Scaling analysis of the routine dgesv
We wish to examine how the execution time grows with $n$, the order of the system of equations to solve. In order to do so,
# Create a large matrix $A$ with random entries. CALL srand(int seed) (seed is an arbitrary integer) will initialize the random sequence, and rand() returns real random numbers uniformly distributed between 0 and 1. Create a matrix of size $(n,n)$ with random entries, and also a right hand side vector $b$. To make the problem better conditioned, add to the diagonal elemens of $A$ some positive quantity (e.g. 1).
# Solve the linear system of equations for increasing values of $n$. Use the time command to obtain the execution time. For example, if a.out is the executable, type instead
{{{
time a.out
}}}
This will give you the execution time. Plot in a log log scale the execution time with $n$ and examine how the time scales with $n$ for large $n$.
Although simple, this method will also include the time involved in the other tasks in your code. For example, it will include the time spent setting up the matrix A. A more accurate timing can be obtained by using the function cpu_time (in fortran):
{{{
CALL cpu_time(t0)
! ... tasks to be timed ...
CALL cpu_time(t1)
WRITE(6,*) 'Elapsed time: ',t1-t0
}}}
In C you would use:
{{{
#include <time.h>
clock_t start, end;
double cpu_time_used;
start = clock();
… /* Do the work. */
end = clock();
cpu_time_used = ((double) (end - start)) / CLOCKS_PER_SEC;
}}}
We will use the LAPACK routine dgesv to solve a partial differential equation.
The heat diffusion equations in one dimension is,
$$ \frac{\partial u(x,t)}{\partial t} = D \frac{\partial^{2} u}{\partial x^{2}}$$
Assume we denote by $u_{i}^{n}=u(x_{i},t_{n})$ the solution at a discrete set of points in some interval of the $x$ axis, $x_{i} = i \Delta x$. $i = 0, \ldots N$, and corresponding to some discrete time $t_{n} = n \Delta t$, with $\Delta x$ and $\Delta t$ the discrete increments in space and time respectively. We are told that the differential equation can be written in these discretized variables as
$$ \frac{1}{\Delta t} \left( u_{i}^{n+1}-u_{i}^{n} \right) = \frac{D}{(\Delta x)^{2}} \left( u_{i+1}^{n+1} - 2 u_{i}^{n+1} + u_{i-1}^{n+1} \right) $$
By rearranging the terms of this equation, an iterative solution can be constructed that relates the solution at some time $t_{n}$, $u_{i}^{n}$, to the solution at the next time step $t_{n+1}$, $u_{i}^{n+1}$. The iteration is,
$$ -a u_{i-1}^{n+1} + (1+2a) u_{i}^{n+1} - a u_{i+1}^{n+1} = u_{i}^{n}$$
with $a = D \Delta t/(\Delta x)^{2}$. This iteration can be written in matrix form $Ax=b$ where $x=\{u_{i}^{n+1} \}$, $b=\{u_{i}^{n} \}$ and the matrix $A$ is $$A = \left( \begin{array}{cccccc} 1+2a & -a & 0 & \ldots& & 0 \\ -a & 1+2a & -a & 0 & \dots & 0 \\ 0 & -a & 1+2a & -a & \ldots & 0 \\ & \ldots & \ldots & \ldots & & \\ 0 & \ldots & 0 & & -a & 1+2a \end{array} \right)$$
As boundary conditions at the end points, we take fixed $u$ for all times: $u_{0} = 0 $ and $u_{N} = 0$.
Write a code that uses the LAPACK routine dgesv to solve the linear system of equations given which amounts to solving the partial differential equation iteratively: Solve for a number of time steps $n = 1, 2, \ldots$ with an initial condition at $t=0$ $u_{k}^{0}=1$ for $k = N/2$ and the remaining $u_{k}^{0} = 0$.
Plot the solution for a few time steps as it evolves in time.
{{justifyright{
''Thursday, September 27, 2018''
}}}
We will use the LAPACK routine dgesv to solve a partial differential equation.
The heat diffusion equations in one dimension is,
$$ \frac{\partial u(x,t)}{\partial t} = D \frac{\partial^{2} u}{\partial x^{2}}$$
Assume we denote by $u_{i}^{n}=u(x_{i},t_{n})$ the solution at a discrete set of points in some interval of the $x$ axis, $x_{i} = i \Delta x$. $i = 0, \ldots N$, and corresponding to some discrete time $t_{n} = n \Delta t$, with $\Delta x$ and $\Delta t$ the discrete increments in space and time respectively. We are told that the differential equation can be written in these discretized variables as
$$ \frac{1}{\Delta t} \left( u_{i}^{n+1}-u_{i}^{n} \right) = \frac{D}{(\Delta x)^{2}} \left( u_{i+1}^{n+1} - 2 u_{i}^{n+1} + u_{i-1}^{n+1} \right) $$
By rearranging the terms of this equation, an iterative solution can be constructed that relates the solution at some time $t_{n}$, $u_{i}^{n}$, to the solution at the next time step $t_{n+1}$, $u_{i}^{n+1}$. The iteration is,
$$ -a u_{i-1}^{n+1} + (1+2a) u_{i}^{n+1} - a u_{i+1}^{n+1} = u_{i}^{n}$$
with $a = D \Delta t/(\Delta x)^{2}$. This iteration can be written in matrix form $Ax=b$ where $x=\{u_{i}^{n+1} \}$, $b=\{u_{i}^{n} \}$ and the matrix $A$ is $$A = \left( \begin{array}{cccccc} 1+2a & -a & 0 & \ldots& & 0 \\ -a & 1+2a & -a & 0 & \dots & 0 \\ 0 & -a & 1+2a & -a & \ldots & 0 \\ & \ldots & \ldots & \ldots & & \\ 0 & \ldots & 0 & & -a & 1+2a \end{array} \right)$$
As boundary conditions at the end points, we take fixed $u$ for all times: $u_{0} = 0 $ and $u_{N} = 0$.
Write a code that uses the LAPACK routine dgesv to solve the linear system of equations given. This amounts to solving the partial differential equation iteratively over $n$. Solve for a number of time steps $n = 1, 2, \ldots$ with an initial condition at $t=0$ $u_{k}^{0}=1$ for $k = N/2$ and the remaining $u_{k}^{0} = 0$.
Plot the solution for a few time steps as it evolves in time.
{{justifyright{
''Thursday, October 1, 2020''
}}}
Both questions below rely on the "trapezoidal rule" of integration, which will be discussed in class. For the purposes of this Lab, it is defined below.
!!!! 1. Integrals over unbounded domains
The trapezoidal rule (and others) which will be discussed in class considers definite integrals over a finite domain $(a,b)$. However, it is possible to either transform the integrands or to analyze convergence with the limits of integration to evaluate integrals over unbounded domains. Propose a suitable plan, and write the associated code to evaluate the following integrals for $s = 0.5, 0.6, 0.7, ..., 3.0$.
## $$\int_{0}^{\infty} \left( x^{3} + sx \right)^{-1/2} dx,$$
## $$\int_{0}^{\infty} \left( x^{2} + 1 \right) ^{-1/2} e^{-sx} dx$$
!!!! 2. Try to evaluate the following integral over an unbounded domain,
$$ I = \int_{\pi}^{\infty} (s + x)^{-1/3} \sin x dx$$ by using the trapezoidal rule, and any of the values of $s$ as above.
@@color(blue):Hint:@@ Compute $I_{n}^{(0)} = \int_{\pi}^{n \pi} (s + x)^{-1/3} \sin x dx$ by using the trapezoidal rule for a range of values of $n$. The resulting sequence $I_{n}^{(0)}$ is alternating and converges very slowly with increasing $n$. To accelerate convergence define $I_{n}^{(1)} = (1/2)(I_{n}^{(0)} + I_{n-1}^{(0)})$. Write a code that iterates this averaging $I_{n}^{(k+1)} = (1/2)(I_{n}^{(k)} + I_{n-1}^{(k)})$ to estimate the value of the integral.
!!!! Trapezoidal rule
Given
$$
\int_{a}^{b} f(x) dx
$$
consider dividing the interval $(a,b)$ in a sequence of equally spaced sub intervals $x_{0} = a, x_{1}, \ldots x_{n} = b$ with $h = x_{1} - x_{0}$ the spacing. Let $f_{i} = f(x_{i})$. Then
$$
\int_{a}^{b} f(x) dx \approx h \left[ \frac{1}{2} f_{0} + f_{1} + \ldots + f_{n-1} + \frac{1}{2} f_{n} \right]
$$
The heat diffusion equations reads, in one dimension,
$$ \frac{\partial u(x,t)}{\partial t} = D \frac{\partial^{2} u}{\partial x^{2}}$$
Assume we denote by $u_{i}^{n}$ the solution at $x = i \Delta x$ corresponding to a time $t = n \Delta t$, with $\Delta x$ and $\Delta t$ the discrete increments in space and time. The differential equation is now discretized as
$$ \frac{1}{\Delta t} \left( u_{i}^{n+1}-u_{i}^{n} \right) = \frac{D}{(\Delta x)^{2}} \left( u_{i+1}^{n+1} - 2 u_{i}^{n+1} + u_{i-1}^{n+1} \right) $$
This iteration can be rewritten as:
$$ -a u_{i-1}^{n+1} + (1+2a) u_{i}^{n+1} - a u_{i+1}^{n+1} = u_{i}^{n}$$
with $a = D \Delta t/(\Delta x)^{2}$. We will also assume boundary conditions at the end points fixed for all times: $u_{0}$ and $u_{N}$.
Write a code that:
* Implements the iteration above in matrix form. Note that the matrix is tridiagonal.
* Implement the recursion relation given in class to solve the resulting tridiagonal system for one time step.
* Solve for a number of time steps in the case in which $u_{0} = u_{N} = 0$, and an initial condition $u_{k}^{0}=1$ for one arbitrary value of $k$ and the remaining $u^{0} = 0$.
!!! Solution of the biharmonic equation in two dimensions
We want to reuse the method that you have developed to solve the Poisson equation to solve the biharmonic equation instead
\begin{equation}
\nabla^{4} V(x,y) + b \nabla^{2} V(x,y) = -\rho(x,y)
\end{equation}
where $b$ is a constant. By introducing the auxiliary function $\nabla^{2} V(x,y) = g(x,y)$, we decompose the original problem into two problems. First, the biharmonic equation can be written as $\nabla^{2}g+bg = -\rho$, which can solved for $g(x,y)$. Then, the definition of $g(x,y)$ is a Poisson equation determining $V(x,y)$.
Use the same discretization scheme and iterative algorithm that you have used in the homework for the Poisson equation (with $V_{ij} = V(x_{i},y_{j})$, $x_{i} = i\Delta x, y_{j} =j \Delta y$):
\begin{equation}
\frac{1}{\Delta x^{2}} \left( V_{i+1,j} - 2 V_{i,j} + V_{i-1,j} + V_{i,j+1} - 2 V_{i,j} + V_{i,j-1} \right) = - \rho_{i,j} \quad i = 1, ..., N; ~ j = 1, ... , M
\end{equation}
# First, modify the algorithm slightly as the equation for $g(x,y)$ is not the Poisson but rather the Helmholtz equation. Use the same iterative method as in the homework to solve this equation.
# Once $g(x,y)$ has been obtained, solve the Poisson equation with the same method to find $V(x,y)$.
{{justifyright{
''Thursday, October 4, 2018''
}}}
!!! Shear stress in turbulent flow
Theory suggests that for turbulent flow in a circular tube of radius $R=1$ the time smoothed velocity field $v_{z}$ along the tube axis is
$$v_{z} = v_{max} \left( 1 - \frac{r}{R} \right)^{\alpha}$$
where $r$ is the radial coordinate inside the tube, and $\alpha$ an unknown exponent. The following measurements have been performed on the velocity as a function of $r$
| $r$ | 0.1 | 0.2 | 0.3 | 0.4 | 0.5 | 0.6 | 0.7 | 0.8 | 0.9 |
| $v_{z}$ | 1.13 | 1.10 | 1.05 | 0.90 | 0.82 | 0.65 | 0.43 | 0.25 | 0.12 |
We are interested in computing from the data the shear stress in the fluid $T$, which in cylindrical coordinates is given by
$$ T = \left( \frac{d^{2} v_{z}}{dr^{2}} + \frac{1}{r} \frac{d v_{z}}{dr} \right)$$
Write a code that computes numerically the shear stress by
# Directly approximating the derivatives by finite differences of the data.
# Interpolate a cubic spline to the data to find $T$ in the interval $r \in (0.1,0.9)$
{{justifyright{
''Thursday, October 8, 2020''
}}}
As an application of the homework problem to find eigenvalues of a matrix, we will use the same method to find the eigenvalues of a differential equation. Imagine that we need to solve the differential equation (with boundary conditions)
$$
\frac{d}{dx} \left[ (1+x^{2}) \frac{dy}{dx} \right] + \lambda y = 0, \quad\quad y(-1)=y(1)=0,
$$
and find the values of $\lambda$ (the eigenvalue) for which this equation with the boundary conditions given has a nontrivial solution ($y = 0$ being a trivial solution).
In order to do it, you can use the following method: Discretize the equation given as,
$$
\frac{d}{dx} \left[ p(x) \frac{dy}{dx} \right] = \frac{1}{h} \left[ p_{n+1/2}\left( \frac{y_{n+1}-y_{n}}{h} \right) - p_{n-1/2} \left( \frac{y_{n}-y_{n-1}}{h} \right) \right] ,
$$
where $p(x) = (1+x^{2})$. Once the differential equation is discretized in this manner, the eigenvalue problem is written in matrix form, and you can find the eigenvalues of the resulting matrix.
Find the smallest eigenvalue.
!!!!1. Continuation of Lab Session 5
* Finish writing the code to solve the diffusion equation as suggested in [[LabSession 5]].
* Integrate que equation forward in time for the following initial conditions:
** A rectangular temperature peak at the center of the domain.
** Two rectangular temperature peaks somewhere in the domain.
** Fixed boundary contitions $u_{0}=2$ and $u_{N} = 1$, with all other $u = 0$ at $t=0$.
** In all three cases, plot the resulting solutions for several times.
!!!!2. Integrals over unbounded domains
Time permitting, try to evaluate the following integral over an unbounded domain,
$$ I = \int_{\pi}^{\infty} (s + x)^{-1/3} \sin x dx$$ by using the trapezoidal rule.
@@color(blue):Hint:@@ Compute $I_{n}^{(0)} = \int_{\pi}^{n \pi} (s + x)^{-1/3} \sin x dx$ by using the trapezoidal rule for a range of values of $n$. The resulting sequence $I_{n}^{(0)}$ is alternating and converges very slowly with increasing $n$. To accelerate convergence define $I_{n}^{(1)} = (1/2)(I_{n}^{(0)} + I_{n-1}^{(0)})$. Write a code that iterates this averaging $I_{n}^{(k+1)} = (1/2)(I_{n}^{(k)} + I_{n-1}^{(k)})$ to estimate the value of the integral.
!!! Integrals over unbounded domains
Evaluate the following integral over an unbounded domain,
$$ I = \int_{\pi}^{\infty} (s + x)^{-1/3} \sin x dx$$ by using the trapezoidal rule.
@@color(blue):Hint:@@ Compute $I_{n}^{(0)} = \int_{\pi}^{n \pi} (s + x)^{-1/3} \sin x dx$ by using the trapezoidal rule for a range of values of $n$. The resulting sequence $I_{n}^{(0)}$ is alternating and converges very slowly with increasing $n$. To accelerate convergence define $I_{n}^{(1)} = (1/2)(I_{n}^{(0)} + I_{n-1}^{(0)})$. Write a code that iterates this averaging $I_{n}^{(k+1)} = (1/2)(I_{n}^{(k)} + I_{n-1}^{(k)})$ to estimate the value of the integral.
!!! Eigenvalue problems for differential equations
Write a code that computes the smallest eigenvalue $\lambda$ for which the problem $$
\frac{d}{dx} \left[ (1+x^{2}) \frac{dy}{dx} \right] + \lambda y = 0, \quad\quad y(-1)=y(1)=0
$$
has a nontrivial solution.
In order to do it, use the following discretization scheme, $$
\frac{d}{dx} \left[ p(x) \frac{dy}{dx} \right] = \frac{1}{h} \left[ p_{n+1/2}\left( \frac{y_{n+1}-y_{n}}{h} \right) - p_{n-1/2} \left( \frac{y_{n}-y_{n-1}}{h} \right) \right] + c_{1}(x) h^{2} + c_{2} h^{4} + \ldots ,
$$
where we have also included the asymptotic error of the discretization. Once the differential equation is discretized, find the smallest eigenvalue of the resulting matrix.
{{justifyright{
''Thursday, October 11, 2018''
}}}
We investigate here a common technique of code development: the use of preprocessors. They "pre process" the source code @@color(red):before@@ compilation. By introducing special commands ("directives") the preprocessor can act on the code prior to compilation, and hence modify the actual code send to the compiler. For that reason, they add a great deal of flexibility in code management and reduce errors by simplifying version maintenance.
!!!Preprocessing (cpp)
As an example, we will write a code that integrates an ordinary differential equation using either forward or backward Euler methods, and we will use a pre-processor to determine which algorithm is used. For C, C++, and FORTRAN code we will use the C preprocessor: cpp (man cpp for help).
cpp processes special statements called directives. They are lines in the source code (in C, C++, and FORTRAN) that begin with #. For example
{{{
# define PI 3.141592
}}}
would literally replace every instance of PI (case sensitive) in a given source code by 3.141592. As an example in FORTRAN, create the code test.f90:
{{{
PROGRAM test
# define PI 3.141592
x = PI
write(6,*) x
END PROGRAM
}}}
In C++ you would use:
{{{
#include <iostream>
#define PI 3.141592
using namespace::std;
int main(){
cout << "pi = " << PI << endl;
return 0;
}
}}}
To run the code through the preprocessor, you would type
{{{
cpp -P test.f90
}}}
The standard ouput will give you the "preprocessed" code. Note how PI has been replaced by its value. To store the preprocessed code in a different file (e.g., test_processed.f90), you would use
{{{
cpp -P -o test_processed.f90 test.f90
}}}
This new file does not contain PI, rather its numerical value. This is a useful feature, for example, in allocating array dimensions without having to edit the source code (even more useful if the modification would involve a large number of modules and subroutines). More about this later.
Another useful feature of defined variables is that they can be used to control code. For example:
{{{
# define LINEAR
....
# if defined(LINEAR)
Execute some instructions
# endif
# if defined(NONLINEAR)
Some other instructions
# endif
}}}
The instructions in the first if block will only appear in the preprocessed code if LINEAR is defined above, but not otherwise. The first line (#define LINEAR) can in effect be used to control whether the LINEAR or NONLINEAR part of the code will be compiled. Try an example yourself, and type cpp -w test.f90 to see the effect.
A third useful directive is
{{{
# include "config.h"
}}}
which includes the file config.h before the preprocessig is done. For example, one could have a file config.h containing:
{{{
# define NSIZE 200
# define LINEAR
}}}
with the code test.f90 including now:
{{{
PROGRAM test
# include "config.h"
REAL(8) x(NSIZE)
# if defined (LINEAR)
do i=1,NSIZE
x(i)=i
END DO
# endif
# if defined(NONLINEAR)
do i=i,NSIZE
x(i)=i**2
END DO
# endif
END PROGRAM
}}}
And the corresponding makefile could read something like
{{{
test: test.f90 config.h
cpp -P -o test_cpp.f90 test.f90
gfortran -o test test_cpp.f90
rm test_cpp.f90
}}}
@@color(red):NOTE:@@ In practice, compilers know how to interpret pre-processor directives without the need to create intermediate files (the files test_cpp.f90 above). In C/C++ the preprocessor is used automatically, and there is nothing for you to do anything at compilation time. In FORTRAN you can either use gfortran -cpp test.f90 or gfortran test.F90. In the first case, the compiler is alerted via -cpp to run the code with the preprocessor. In the second, the convention is used that if the extension has a capital F (test.F90) not a lower case (test.f90), the compiler will use the preprocessor first.
Use this method to write a code that solves the ordinary differential equation
$$ \frac{dy}{dt} = \frac{\sin y}{\sqrt{y}} $$
with initial condition at $t=0$ $y=1$ up to $t=1$. Write a code that can use either the forward or backward Euler method and use the pre-processor to choose which method to use. For example, create a config.h file that contains # define FORWARD or # define BACKWARD and the corresponding sections of the main code (e.g. euler.f90 or euler.c) to use one algorithm or the other depending on the contents of the file config.h
{{justifyright{
''Thursday, October 15, 2020''
}}}
!!! Study of planetary motion
According to Newton's Law of gravitation, the force on a planet of mass $m$ orbiting the sun of mass $M$ is
$$
F = - \frac{GMm}{r^{2}},
$$
where the force is directed along the radial direction, $\vec{r}$ is the vector joining the Sun and the planet, and $r = \| \vec{r} \|$. $G$ is the universal gravitation constant. Given Newton's Second Law, and denoting $\vec{r} = (x,y)$, the two coordinates of the planet, the differential equations governing its motion are,
$$
\frac{d^{2}x}{dt^{2}} = - GM \frac{x}{r^{3}}, \quad\quad \frac{d^{2}y}{dt^{2}} = - GM \frac{y}{r^{3}}
$$
Assume units in which $GM = 1$, and use the initial conditions,
$$
x(0) = 0.5, \quad y(0) = 0, \quad v_{x}(0) = 0, \quad v_{y}(0) = 1.63
$$
# Integrate these equations over time to complete one orbit. You will need to take the time step small enough; otherwise the orbit will not close onto itself.
# Experiment with the initial conditions until you find a circular orbit.
# Once the algorithm and the time step lead to good accuracy, progressively increase the initial velocity until the orbits cease being elliptic and become open (hyperbolic - the planet flies away).
# Using the same initial conditions as above, examine the effect of using a $1/r^{4}$ power in the law of gravitation (instead of $1/r^{2}$). You should find that the ellipse now rotates (precession). In fact, you can verify that any slight variation of the exponent away from 2 leads to precession.
!!!!Boundary value problem
We wish to solve the following second order equation in one dimension
$$y^{\prime\prime} = -\frac{4}{x} y^{\prime} + \frac{2}{x^{2}} y - \frac{2}{x^{2}} \ln x$$
in the interval $1 \le x \le 2$, and with boundary conditions $y(1) = -1/2$, $y(2) = \ln 2$.
# Write the second order equation as a system of two first order equations, each with one boundary condition.
# Use a method of your choice to discretize and solve the system of two equations as if it were an initial value problem.
# We will use a "shooting method" to solve the boundary value problem. Set up the solution method by assuming one missing boundary condition at $x = 1$ and integrate the system up to $x=2$.
# Set up an interative method so that changes in the assumed boundary condition at $x=1$ lead, after convergence, to the correct boundary condition at $x=2$, and hence to the solution of the boundary value problem.
A simple minded solution. The code bvp.f90:
{{{
PROGRAM bvp
p0_old=3.0
dp=0.01
CALL euler(0.01,-0.5,p0_old,yn1,pn1)
p0_plus=p0_old+dp
CALL euler(0.01,-0.5,p0_plus,yn2,pn2)
p0_new=p0_old-(yn1-log(2.0))*dp/(yn2-yn1)
CALL euler(0.01,-0.5,p0_new,yn3,pn3)
WRITE(6,*) p0_old,yn1-log(2.0),p0_new,yn3-log(2.0)
END PROGRAM
}}}
whis uses the subroutine euler.f90
{{{
SUBROUTINE euler(dx,y0,p0,yn,pn)
open(1,file='output.dat')
yn=y0
pn=p0
niter=(2.0-1.0)/dx
do i=1,niter
ysave=yn
x=1.0+i*dx
yn=yn+dx*pn
pn=pn+dx*(-4.0*pn/x+2.0*ysave/x**2-2.0*log(x)/x**2)
write(1,*) x,yn,pn
end do
CLOSE(1)
RETURN
END SUBROUTINE
}}}
and the makefile to put it all together
{{{
bvp: bvp.f90 euler.o
gfortran -o bvp bvp.f90 euler.o
euler.o: euler.f90
gfortran -c euler.f90
}}}
Numerical integration of a predator-prey model. Consider two species of fish that coexist in the same lake. Let $x$ be the
!!! Numerical integration of a predator-prey model.
Consider two species of fish that coexist in the same lake. Let $x$ be number of small fish of type A which only eat grass from the bottom, and $y$ the number of large fish of type B which only eat small fish of type A, but no grass. Alone, the population of A might be expected to grow exponentially, whereas B alone (no A to eat) would die exponentially. However, the two species interact: whenever a fish B encounters a fish A, it will eat it. This interaction contributes to the increase in the population of B and a decrease in A, both rates being proportional to the quantity $xy$. A general mathematical model of this system is $$ \frac{dx}{dt} = \epsilon_{A} x - \gamma_{A} xy $$
$$ \frac{dy}{dt} = - \epsilon_{B} y + \gamma_{B} xy $$ where $\epsilon_{A}$ is the growth rate of A alone, $\epsilon_{B}$ the extinction of B alone, and $\gamma_{A}$ and $\gamma_{B}$ the interaction coefficients. All four are assumed to be constant. Write a code that integrates this system of differential equations with $\epsilon_{A} = \epsilon_{B} = \gamma_{A} = \gamma_{B} =1$, and initial conditions at $t=0$ $x=0.1$ and $y=0.1$.
## Integrate the system up to $t=10$.
## Plot $y$ versus $x$ and provide an interpretation of the resulting graph.
{{justifyright{
''Thursday, October 18, 2018''
}}}
!!! Boundary value problem
We wish to solve the following second order equation in one dimension
$$y^{\prime\prime} = -\frac{4}{x} y^{\prime} + \frac{2}{x^{2}} y - \frac{2}{x^{2}} \ln x$$
in the interval $1 \le x \le 2$, and with boundary conditions $y(1) = -1/2$, $y(2) = \ln 2$.
# Write the second order equation as a system of two first order equations, each with one boundary condition.
# Use a method of your choice to discretize and solve the system of two equations as if it were an initial value problem.
# We will use a "shooting method" to solve the boundary value problem. Set up the solution method by assuming one missing boundary condition at $x = 1$ and integrate the system up to $x=2$.
# Set up an interative method so that changes in the assumed boundary condition at $x=1$ lead, after convergence, to the correct boundary condition at $x=2$, and hence to the solution of the boundary value problem.
{{justifyright{
''Thursday, October 22, 2020''
}}}
!!! Classical scattering problem
We wish to study the two dimensional motion of a particle that scatters off the following potential
$$
V(x,y) = x^{2}y^{2} e^{-(x^{2}+y^{2})}.
$$
Begin by plotting the potential to get an idea of the behavior. You will see four peaks in potential symmetric around $(0,0)$. A scattering experiment involves a particle starting far away, and moving towards the origin. An experiment would record the location and orientation of the scattered particle. In order to solve this problem, we write Newton's equations of motion,
$$
m \frac{d^{2}x}{dt^{2}} = - 2y^{2}x(1-x^{2}) e^{-(x^{2}+y^{2})},
$$
$$
m \frac{d^{2}y}{dt^{2}} = - 2x^{2}y(1-y^{2}) e^{-(x^{2}+y^{2})}.
$$
Good parameters for the numerical solution are $m=0.5$, and initial conditions $v_{x}(0) = 0, v_{y}(0) = 0.5$. The initial position of the particle $y(0)$ is large negative $y$. Make it large enough so that the exponential terms in the potential are negligible. Explore a range of initial $b = x(0) = (-1,1)$ (this value $b$ is called the "impact parameter"). With this condition, the particle will approach the origin and scatter off the four peaks.
# Plot a number of trajectories $(x(t),y(t))$ for different values of $b$ that show usual and unusual behavior (e.g., backward scattering, multiple scattering, etc.).
# For the trajectories obtained, study "phase space plots", that is, plots with axes $(x(t), dx(t)/dt)$ and $(y(t), dy(t)/dt)$. How do these trajectories differ from "bound states" (those in which the particle remains trapped near $(x,y)=(0,0)$).
# A "scattering angle" is defined as $\theta =$atan$(vy/vx)$ far from the origin (use the function {{{theta=atan2(vx,vy)}}} on the computer). Plot the scattering angle as a function of $b$ and note any anomalies.
We investigate here a common technique of code development: the use of preprocessors. They "pre process" the source code @@color(red):before@@ compilation. By introducing special commands ("directives") the preprocessor can act on the code prior to compilation, and hence modify the actual code send to the compiler. For that reason, they add a great deal of flexibility in code management and reduce errors by simplifying version maintenance.
!!!Preprocessing (cpp)
As an example, we will write a code that integrates an ordinary differential equation using either forward or backward Euler methods, and we will use a pre-processor to determine which algorithm is used. For both C and FORTRAN code we will use the C preprocessor: cpp (man cpp for help).
cpp processes special statements called directives. They are lines in the source code (both C and FORTRAN) that begin with #. For example
{{{
# define PI 3.141592
}}}
would literally replace every instance of PI (case sensitive) in a given source code by 3.141592. As an example, create the code test.f90:
{{{
PROGRAM test
# define PI 3.141592
x = PI
write(6,*) x
END PROGRAM
}}}
To run the code through the preprocessor, you would type
{{{
cpp -w test.f90
}}}
The standard ouput will give you the "preprocessed" code. Note how PI has been replaced by its value. To store the preprocessed code in a different file (e.g., test_processed.f90), you would use
{{{
cpp -w -o test_processed.f90 test.f90
}}}
This new file does not contain PI, rather its numerical value. This is a useful feature, for example, in allocating array dimensions without having to edit the source code (even more useful if the modification would involve a large number of modules and subroutines). More about this later.
Another useful feature of defined variables is that they can be used to control code. For example:
{{{
# define LINEAR
....
# if defined(LINEAR)
Execute some instructions
# endif
# if defined(NONLINEAR)
Some other instructions
# endif
}}}
The instructions in the first if block will only appear in the preprocessed code if LINEAR is defined above, but not otherwise. The first line (#define LINEAR) can in effect be used to control whether the LINEAR or NONLINEAR part of the code will be compiled. Try an example yourself, and type cpp -w test.f90 to see the effect.
A third useful directive is
{{{
# include "config.h"
}}}
which includes the file config.h before the preprocessig is done. For example, one could have a file config.h containing:
{{{
# define NSIZE 200
# define LINEAR
}}}
with the code test.f90 including now:
{{{
PROGRAM test
# include "config.h"
REAL(8) x(NSIZE)
# if defined (LINEAR)
do i=1,NSIZE
x(i)=i
END DO
# endif
# if defined(NONLINEAR)
do i=i,NSIZE
x(i)=i**2
END DO
# endif
END PROGRAM
}}}
And the corresponding makefile could read something like
{{{
test: test.f90 config.h
cpp -w -o test_cpp.f90 test.f90
gfortran -o test test_cpp.f90
rm test_cpp.f90
}}}
Use this method to write a code that solves the ordinary differential equation
$$ \frac{dy}{dt} = \frac{\sin y}{\sqrt{y}} $$
with initial condition at $t=0$ $y=1$ up to $t=1$. Write a code that can use either the forward or backward Euler method and use the pre-processor to choose which method to use. For example, create a config.h file that contains # define FORWARD or # define BACKWARD and the corresponding sections of the main code (e.g. euler.f90 or euler.c) to use one algorithm or the other depending on the contents of the file config.h
!!!!Boundary value problem
We wish to solve the following second order equation in one dimension
$$y^{\prime\prime} = -\frac{4}{x} y^{\prime} + \frac{2}{x^{2}} y - \frac{2}{x^{2}} \ln x$$
in the interval $1 \le x \le 2$, and with boundary conditions $y(1) = -1/2$, $y(2) = \ln 2$.
# Write the second order equation as a system of two first order equations, each with one boundary condition.
# Use a method of your choice to discretize and solve the system of two equations as if it were an initial value problem.
# We will use a "shooting method" to solve the boundary value problem. Set up the solution method by assuming one missing boundary condition at $x = 1$ and integrate the system up to $x=2$.
# Set up an interative method so that changes in the assumed boundary condition at $x=1$ lead, after convergence, to the correct boundary condition at $x=2$, and hence to the solution of the boundary value problem.
A simple minded solution. The code bvp.f90:
{{{
PROGRAM bvp
p0_old=3.0
dp=0.01
CALL euler(0.01,-0.5,p0_old,yn1,pn1)
p0_plus=p0_old+dp
CALL euler(0.01,-0.5,p0_plus,yn2,pn2)
p0_new=p0_old-(yn1-log(2.0))*dp/(yn2-yn1)
CALL euler(0.01,-0.5,p0_new,yn3,pn3)
WRITE(6,*) p0_old,yn1-log(2.0),p0_new,yn3-log(2.0)
END PROGRAM
}}}
whis uses the subroutine euler.f90
{{{
SUBROUTINE euler(dx,y0,p0,yn,pn)
open(1,file='output.dat')
yn=y0
pn=p0
niter=(2.0-1.0)/dx
do i=1,niter
ysave=yn
x=1.0+i*dx
yn=yn+dx*pn
pn=pn+dx*(-4.0*pn/x+2.0*ysave/x**2-2.0*log(x)/x**2)
write(1,*) x,yn,pn
end do
CLOSE(1)
RETURN
END SUBROUTINE
}}}
and the makefile to put it all together
{{{
bvp: bvp.f90 euler.o
gfortran -o bvp bvp.f90 euler.o
euler.o: euler.f90
gfortran -c euler.f90
}}}
{{justifyright{
''Thursday, October 25, 2018''
}}}
!!! ~Runge-Kutta integration of a system of ordinary differential equations
We wish to calculate the vertical density distribution of a self-gravitating plane distribution of gas in isothermal, hydrostatic equilibrium by using the ~Runge-Kutta solver of order $4^{th}$. The model is defined by by the following two, coupled ordinary differential equations,
$$
\frac{d \rho(z)}{dz} = \frac{g}{c^{2}} \rho ,
$$
$$
\frac{d g(z)}{dz} = - 4 \pi G \rho ,
$$
where $\rho(z)$ is the mass density and $g(z)$ the self-consistently determined intensity of the gravitational field. $c^{2}$ and $G$ are known constants. For convenience, you can set $c^{2} = 2 \pi G = 1$.
# Write an algorithm using the ~Runge-Kutta method of order 4 to solve this system of equations with boundary conditions $\rho(z=0) = \rho_{0} = 1$ and $g(z=0) =0$ up to $z_{max} = 10$. A discretization step $\Delta z = 0.4$ should be sufficient.
# For reference, compare your solution to the known analytic solution given by, $$
\rho(z) = \rho_{0} \; {\rm sech}^{2}(z/H) \quad g(z) = -4 \pi G \Sigma \; {\rm tanh}(z/H)
$$
where $H = c^{2}/(2 \pi G \Sigma)$ and $\Sigma = \int_{0}^{\infty} \rho(z) dz$. You can estimate $\Sigma$ by using the trapezoidal rule on the numerically determined density profile.
Refresh your git repository with updated tridiagonal routines. Go to the directory in which you have the respository and type
{{{
git pull
}}}
!!! Alternating Direction Implicit solution of Poisson's equation in two dimensions
As we have seen in class, a simple discretization of Poisson's equations on a two dimensional square grid of spacing $\Delta x$ is
$$u_{i+1 j} + u_{i-1 j} + u_{i j-1} + u_{i j+1} - 4 u_{ij} = - \Delta x^{2} \rho_{ij}$$
In order to develop and ADI algorithm, we split it as follows
$$\left( u_{i+1 j} + u_{i-1 j} - 2 u_{ij} \right) + \left( u_{i j-1} + u_{i j+1} - 2 u_{ij} \right) = - \Delta x^{2} \rho_{ij}$$
In the first parenthesis the index $j$ is constant, whereas in the second parenthesis the index $i$ is constant.
!!!!First half step
We define the first half step of the ADI method assuming implicit updating in $i$ only:
$$u_{i+1 j}^{n+1/2} + u_{i-1 j}^{n+1/2} - 2 u_{ij}^{n+1/2} = - \left( u_{i j-1}^{n} + u_{i j+1}^{n} - 2 u_{ij}^{n} \right) - \Delta x^{2} \rho_{ij}$$
@@color(maroon):If you now think that $j$ is fixed@@, then for each $j$ this equations is just a standard tri-diagonal system - or in effect $n$ tridiagonal systems, one for each $j$. The right hand side depends on $\rho_{ij}$ but also on the previous solution $u_{ij}$. Note that the right hand side involves different values of $j$, but since they are at an earlier time, they are known.
We also need to specify boundary conditions. For now, let us take $u = 0$ on all boundaries. Then, in the notation given in class, the elements of the tridiagonal matrix are:
\begin{eqnarray}
b(1) = b(n) = 1 & \quad & b(2) ... b(n-1) = -2 \\
a(n) = 0 & \quad & a(2) ... a(n-1) = 1\\
c(1) = 0 & \quad & c(2) ... c(n-1) = 1
\end{eqnarray}
Write a loop over $j$ so that for each $j$ the tridiagonal system above is solved, leading to $u_{ij}^{n+1/2}$. Use the routine provided: tridag.f or tridag.c.
!!!!Second half step
We now reverse the implicit direction for the second half step:
$$u_{i j+1}^{n+1} + u_{i j-1}^{n+1} - 2 u_{ij}^{n+1} = - \left( u_{i+1 j}^{n+1/2} + u_{i-1 j}^{n+1/2} - 2 u_{ij}^{n+1/2} \right) - \Delta x^{2} \rho_{ij}$$
@@color(maroon):If you now think of $i$ being fixed instead@@, for each $i$ we have a tridiagonal system on $j$. So in effect, we have again $n$ tridiagonal systems. Furthermore, note that the form of the tridiagonal system is identically the same as above.
Write a second loop to solve all $n$ tridiagonal systems.
This completes one iteration in the solution to Poisson's equation for any given $\rho_{ij}$ and boundary conditions. The process can be iterated until convergence. Sample code:
{{{
PROGRAM poisson
PARAMETER (n=100)
REAL(8) :: u(n,n),un(n,n),rho(n,n),r(n,n)
REAL(8) :: a(n),b(n),c(n)
dx=0.1
do i=1,n
do j=1,n
rho(i,j)=-dx**2*sqrt(FLOAT(i)**2+FLOAT(j)**2)
u(i,j)=0.0
un(i,j)=0.0
end do
end do
a(n)=0.0
do i=1,n-1
a(i)=1.0
end do
b(1)=1.0
b(n)=1.0
do i=2,n-1
b(i)=-2.0
end do
c(1)=0.0
do i=2,n-1
c(i)=1.0
end do
DO j=2,n-1
DO i=1,n
r(i,j)=rho(i,j)-(u(i,j+1)+u(i,j-1)-2.0*u(i,j))
END DO
call tridag(a,b,c,r(:,j),un(:,j),n)
END DO
DO i=2,n-1
DO j=1,n
r(i,j)=rho(i,j)-(un(i+1,j)+un(i-1,j)-2.0*un(i,j))
END DO
call tridag(a,b,c,r(i,:),u(i,:),n)
END DO
WRITE(6,*) 'One iteration completed'
END PROGRAM
}}}
!!! Alternating Direction Implicit (ADI) method for convection-diffusion equation
The temperature distribution in a porous rock in the presence of ground water flow directed along the $x$ axis is described by the equation
$$ v \frac{\partial T}{\partial x} = \kappa \left( \frac{\partial^{2} T}{\partial x^{2}} + \frac{\partial^{2} T}{\partial y^{2}} \right) $$
where $v$ is the constant water velocity and $\kappa$ is the thermal diffusivity coefficient.
Each iteration of the ADI method involves two consecutive steps, each being semi-implicit in alternating directions. In this case, considering central differences for the spatial derivatives, the two steps are $$
(1 -\beta) T_{i+1,j}^{n+1/2} + 2 \beta T_{i,j}^{n+1/2} - (1+\beta) T_{i-1,j}^{n+1/2} = \beta \left( T_{i,j+1}^{n} -2 T_{i,j}^{n} + T_{i,j-1}^{n} \right)
$$
and $$
\beta \left( T_{i,j+1}^{n+1} -2 T_{i,j}^{n+1} + T_{i,j-1}^{n+1} \right) = (1-\beta) T_{i+1,j}^{n+1/2} + 2 \beta T_{i,j}^{n+1/2} - (1+\beta) T_{i-1,j}^{n+1/2}
$$
with $\beta = 2 \kappa/v \Delta x$ ($\Delta y = \Delta x$ here). The first step runs over all values of $j = 1, \ldots n$, and the second step runs over all values of $i = 1, \ldots, n$. Both are a tri-diagonal system of linear equations.
Sketch the solution to the two systems of linear equations defined by these two iterations by using the routine tridag (FORTRAN or C) uploaded to the class repository. As you will see, the routine takes as input three arrays: $b_{k}, k=1, \ldots N = n \times n \;$ contains the diagonal elements in the matrix representing the system, $a_{k}, k =2, \ldots, N$ the line immediately below the diagonal, and $c_{k}, k=1, \ldots N-1$ the line immediately above the diagonal. The right hand side is the array $r_{k}, k=1, \ldots N$.
In order to do so, map the two ADI steps given above to the corresponding arrays $a,b,c, r$. Then, given the output $u_{k}, k=1, \ldots N$ of tridag, map it back to the temperature $T_{ij}$.
{{justifyright{
''Thursday, November 1, 2018''
}}}
!!! Soliton solution
We wish to find nonlinear solitary wave solutions ("solitons") of the following equation $$
\frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} + \delta^{2} \frac{\partial^{3} u}{\partial x^{3}} = 0.
$$
where we choose $\delta = 0.022$ and take as initial condition $u(x,0) = \cos (\pi x)$, $0 \le x \le 2$. We will also impose boundary conditions that $u$ is periodic in $[0,2]$ for all time.
Develop a code that uses the ~Lax-Wendroff method to integrate the equation in time, and examine the emergence of the soliton solution. Note that you first have to write the equation given in conservative form $\partial_{t} u + \partial_{x}J(u) = 0$. Integrate the equation up to a maximum time $t=3.6/\pi$.
All lectures have been recorded. You can find the recordings up to October 12, 2020 [[here|https://canvas.umn.edu/courses/195750/files/folder/Recorded%20Lectures]]. From October 13, 2020, you can find the recordings [[here|https://canvas.umn.edu/courses/195750/external_tools/23]].
* [[September 2020|Schedule_SEP2020]].
* [[October 2020|Schedule_OCT2020]].
/***
|''Name:''|LoadRemoteFileThroughProxy (previous LoadRemoteFileHijack)|
|''Description:''|When the TiddlyWiki file is located on the web (view over http) the content of [[SiteProxy]] tiddler is added in front of the file url. If [[SiteProxy]] does not exist "/proxy/" is added. |
|''Version:''|1.1.0|
|''Date:''|mar 17, 2007|
|''Source:''|http://tiddlywiki.bidix.info/#LoadRemoteFileHijack|
|''Author:''|BidiX (BidiX (at) bidix (dot) info)|
|''License:''|[[BSD open source license|http://tiddlywiki.bidix.info/#%5B%5BBSD%20open%20source%20license%5D%5D ]]|
|''~CoreVersion:''|2.2.0|
***/
//{{{
version.extensions.LoadRemoteFileThroughProxy = {
major: 1, minor: 1, revision: 0,
date: new Date("mar 17, 2007"),
source: "http://tiddlywiki.bidix.info/#LoadRemoteFileThroughProxy"};
if (!window.bidix) window.bidix = {}; // bidix namespace
if (!bidix.core) bidix.core = {};
bidix.core.loadRemoteFile = loadRemoteFile;
loadRemoteFile = function(url,callback,params)
{
if ((document.location.toString().substr(0,4) == "http") && (url.substr(0,4) == "http")){
url = store.getTiddlerText("SiteProxy", "/proxy/") + url;
}
return bidix.core.loadRemoteFile(url,callback,params);
}
//}}}
[[Home]]
[[Syllabus]]
[[Lecture schedule]]
[[Homeworks, Labs and Exams]]
[[Supplementary Reading]]
[[Tag Cloud]]
[[WelcomeToTiddlyspot]] [[GettingStarted]]
{{justifyright{
''Due Wednesday, October 22, 2014''
}}}
!!!!Shear stress in turbulent flow
Theory suggests that for turbulent flow in a circular tube of radius $R=1$ the time smoothed velocity field $v_{z}$ along the tube axis is
$$v_{z} = v_{max} \left( 1 - \frac{r}{R} \right)^{\alpha}$$
where $r$ is the radial coordinate inside the tube, and $\alpha$ an unknown exponent. The following measurements have been performed on the velocity as a function of $r$
| $r$ | 0.1 | 0.2 | 0.3 | 0.4 | 0.5 | 0.6 | 0.7 | 0.8 | 0.9 |
| $v_{z}$ | 1.13 | 1.10 | 1.05 | 0.90 | 0.82 | 0.65 | 0.43 | 0.25 | 0.12 |
We are interested in computing from the data the shear stress in the fluid $T$, which in cylindrical coordinates is given by
$$ T = \left( \frac{d^{2} v_{z}}{dr^{2}} + \frac{1}{r} \frac{d v_{z}}{dr} \right)$$
Write a code that computes numerically the shear stress by
# Directly approximating the derivatives by finite differences of the data.
# Interpolate a cubic spline to the data to find $T$ in the interval $r \in (0.1,0.9)$
!!!!Differential equation solution
We wish to solve the equation
$$ \frac{dy}{dt} = \frac{\sin y}{\sqrt{y}} $$
with initial condition at $t=0$ $y=1$, with a predictor-corrector method. For the explicit part we will use the forward Euler method, whereas the implicit part of the algorithm uses the trapezoidal rule
$$ y_{i+1} = y_{i} + \frac{\Delta t}{2} \left[ f(t_{i},y_{i}) + f(t_{i+1},y_{i+1}) \right] $$
where $f(t,y)$ is the right hand side of the differential equation.
# Take $\Delta t = 0.1$ and write a code that calculates the solution up to $t=1$. For the implicit part, set up a fixed point iteration at each step to find the solution up to six decimal places.
# Derive the convergence condition for the fixed point iteration used, and give an upper bound for the time step $\Delta t$ that is necessary for the iteration to converge.
[[pdf version|http://homepages.spa.umn.edu/~vinals/tspot_files/phys4041/2014/PHYS4041_MidTerm.pdf]]
{{justifyright{
''Due Friday, October 21, 2016, 5:00 p.m.''
}}}
We need to find the solution of the following integral equation in the interval $(0,1)$,$$
\int_{0}^{1} dx e^{-(x-y)^{2}/\epsilon^{2}} \frac{df(x)}{dx} = 2 \epsilon \sqrt{\pi} ( y - 3 y^{2} + 2 y^{3}),
$$
where the function $f(x)$ needs to satisfy the boundary conditions $f(0)=f(1)=0$. Also, $\epsilon = 0.1$.
Discretize $x_{j}= j h, j = 0, \ldots n$ for some $h$, so that $x_{0} = 0$ and $x_{n} = 1$, and $y_{i} = i h, i = 0, \ldots n$. Consistent with this discretization, one defines $f_{j} = f(x_{j})$. The boundary conditions are then $f_{0} = f_{n} = 0$.
# Write the discretized version of the integral equation by using the trapezoidal rule for the integral and central differences for the derivative at the inner points, the forward difference at $x=0$ and the backward difference at $x=1$.
# Use the discretized version to write the code that will find the solution $f(x)$ in $(0,1)$.
{{justifyright{
''Due Friday, October 19, 2018, 5:00 p.m.''
}}}
@@color(red):''NOTE:''@@ Just to re iterate that you are expected to work //independently// when solving the mid term exam. You can submit the exam electronically as you have done for the homework.
!!! Question 1. Numerical integration over an unbounded domain (50 points).
The integration methods discussed in class consider definite integrals over a finite domain $(a,b)$. However, it is possible to either transform the integrands or to analyze convergence with the limits of integration to evaluate integrals over unbounded domains. Propose a suitable plan, and write the associate code to evaluate the following integrals for $s = 0.5, 0.6, 0.7, ..., 3.0$.
## $$\int_{0}^{\infty} \left( x^{3} + sx \right)^{-1/2} dx,$$
## $$\int_{0}^{\infty} \left( x^{2} + 1 \right) ^{-1/2} e^{-sx} dx$$
!!! Question 2. Two predators, one prey model (50 points).
Consider two predator species of numbers $P_{1}$ and $P_{2}$ that prey on the same species of number $p$. A possible model of evolution in time $t$ given the interaction between the three species is,
$$
\frac{dp}{dt} = ap \left( 1 - \frac{p}{K} \right) - \left( b_{1}P_{1} + b_{2}P_{2} \right) p,
$$
$$
\frac{dP_{1}}{dt} = \epsilon_{1} b_{1} pP_{1} - m_{1} P_{1}, \quad\quad \frac{dP_{2}}{dt} = \epsilon_{2} b_{2} pP_{2} - m_{2} P_{2}.
$$
Consider the following values of the model parameters: $a = 0.2, K = 1.7, b_{1} = 0.1, b_{2} = 0.2, m_{1} = m_{2} = 0.1, \epsilon_{1}=1.0, \epsilon_{2} = 2.0$.
# Integrate the model forward in time to determine the populations as a function of time, given the initial conditions $p(0) = P_{2}(0) = 1.7$ and $P_{1}(0) = 1.0$.
# To understand the time dependence of the model, it may help to not only plot the functions $p(t), P_{1}(t)$ and $P_{2}(t)$ as a function of time, but also against each other.
# Based on the results obtained for the evolution, can two predators that share the same prey coexist ? You may want to change the parameters of the second predator and examine the changes induced in the evolution.
# Search for chaotic (random) behavior by varying the growth rates $a, \epsilon_{1}$ and $\epsilon_{2}$.
[[pdf version|http://homepages.spa.umn.edu/~vinals/tspot_files/phys4041/2018/PHYS4041_MidTerm.pdf]]
/***
|''Name:''|PasswordOptionPlugin|
|''Description:''|Extends TiddlyWiki options with non encrypted password option.|
|''Version:''|1.0.2|
|''Date:''|Apr 19, 2007|
|''Source:''|http://tiddlywiki.bidix.info/#PasswordOptionPlugin|
|''Author:''|BidiX (BidiX (at) bidix (dot) info)|
|''License:''|[[BSD open source license|http://tiddlywiki.bidix.info/#%5B%5BBSD%20open%20source%20license%5D%5D ]]|
|''~CoreVersion:''|2.2.0 (Beta 5)|
***/
//{{{
version.extensions.PasswordOptionPlugin = {
major: 1, minor: 0, revision: 2,
date: new Date("Apr 19, 2007"),
source: 'http://tiddlywiki.bidix.info/#PasswordOptionPlugin',
author: 'BidiX (BidiX (at) bidix (dot) info',
license: '[[BSD open source license|http://tiddlywiki.bidix.info/#%5B%5BBSD%20open%20source%20license%5D%5D]]',
coreVersion: '2.2.0 (Beta 5)'
};
config.macros.option.passwordCheckboxLabel = "Save this password on this computer";
config.macros.option.passwordInputType = "password"; // password | text
setStylesheet(".pasOptionInput {width: 11em;}\n","passwordInputTypeStyle");
merge(config.macros.option.types, {
'pas': {
elementType: "input",
valueField: "value",
eventName: "onkeyup",
className: "pasOptionInput",
typeValue: config.macros.option.passwordInputType,
create: function(place,type,opt,className,desc) {
// password field
config.macros.option.genericCreate(place,'pas',opt,className,desc);
// checkbox linked with this password "save this password on this computer"
config.macros.option.genericCreate(place,'chk','chk'+opt,className,desc);
// text savePasswordCheckboxLabel
place.appendChild(document.createTextNode(config.macros.option.passwordCheckboxLabel));
},
onChange: config.macros.option.genericOnChange
}
});
merge(config.optionHandlers['chk'], {
get: function(name) {
// is there an option linked with this chk ?
var opt = name.substr(3);
if (config.options[opt])
saveOptionCookie(opt);
return config.options[name] ? "true" : "false";
}
});
merge(config.optionHandlers, {
'pas': {
get: function(name) {
if (config.options["chk"+name]) {
return encodeCookie(config.options[name].toString());
} else {
return "";
}
},
set: function(name,value) {config.options[name] = decodeCookie(value);}
}
});
// need to reload options to load passwordOptions
loadOptionsCookie();
/*
if (!config.options['pasPassword'])
config.options['pasPassword'] = '';
merge(config.optionsDesc,{
pasPassword: "Test password"
});
*/
//}}}
/***
|Name|Plugin: Scientific Notation|
|Created by|BobMcElrath|
|Email|my first name at my last name dot org|
|Location|http://bob.mcelrath.org/tiddlyjsmath-2.0.3.html|
|Version|1.0|
|Requires|[[TiddlyWiki|http://www.tiddlywiki.com]] ≥ 2.0.3, [[jsMath|http://www.math.union.edu/~dpvc/jsMath/]] ≥ 3.0, [[Plugin: jsMath]]|
!Description
This plugin will render numbers expressed in scientific notation, such as {{{3.5483e12}}} using the jsMath plugin to display it in an intuitive way such as 3.5483e12. You may customize the number of significant figures displayed, as well as "normalize" numbers so that {{{47392.387e9}}} is displayed as 47392.387e9.
!Installation
Install the Requirements, above, add this tiddler to your tiddlywiki, and give it the {{{systemConfig}}} tag.
!History
* 1-Feb-06, version 1.0, Initial release
!Code
***/
//{{{
config.formatters.push({
name: "scientificNotation",
match: "\\b[0-9]+\\.[0-9]+[eE][+-]?[0-9]+\\b",
element: "span",
className: "math",
normalize: true, // set to 'true' to convert numbers to X.XXX \times 10^{y}
sigfigs: 3, // with this many digits in the mantissa
handler: function(w) {
var snRegExp = new RegExp("\\b([0-9]+(?:\\.[0-9]+)?)[eE]([-0-9+]+)\\b");
var mymatch = snRegExp.exec(w.matchText);
var mantissa = mymatch[1];
var exponent = parseInt(mymatch[2]);
// normalize the number.
if(this.normalize) {
mantissa = parseFloat(mantissa);
while(mantissa > 10.0) {
mantissa = mantissa / 10.0;
exponent++;
}
while(mantissa < 1.0) {
mantissa = mantissa * 10.0;
exponent--;
}
var sigfigsleft = this.sigfigs;
mantissa = parseInt(mantissa) + "." + (Math.round(Math.pow(10,this.sigfigs-1)*mantissa)+"").substr(1,this.sigfigs-1);
}
var e = document.createElement(this.element);
e.className = this.className;
if(exponent == 0) {
e.appendChild(document.createTextNode(mantissa));
} else {
e.appendChild(document.createTextNode(mantissa + "\\times 10^{" + exponent + "}"));
}
w.output.appendChild(e);
}
});
//}}}
/***
|Name|Plugin: arXiv Links|
|Created by|BobMcElrath|
|Email|my first name at my last name dot org|
|Location|http://bob.mcelrath.org/tiddlyjsmath-2.0.3.html|
|Version|1.0|
|Requires|[[TiddlyWiki|http://www.tiddlywiki.com]] ≥ 2.0.3|
!Description
This formatting plugin will render links to the [[arXiv|http://www.arxiv.org]] preprint system. If you type a paper reference such as hep-ph/0509024, it will be rendered as an external link to the abstract of that paper.
!Installation
Add this tiddler to your tiddlywiki, and give it the {{{systemConfig}}} tag.
!History
* 1-Feb-06, version 1.0, Initial release
!Code
***/
//{{{
config.formatters.push({
name: "arXivLinks",
match: "\\b(?:astro-ph|cond-mat|hep-ph|hep-th|hep-lat|gr-qc|nucl-ex|nucl-th|quant-ph|(?:cs|math|nlin|physics|q-bio)(?:\\.[A-Z]{2})?)/[0-9]{7}\\b",
element: "a",
handler: function(w) {
var e = createExternalLink(w.output, "http://arxiv.org/abs/"+w.matchText);
e.target = "_blank"; // open in new window
w.outputText(e,w.matchStart,w.nextMatch);
}
});
//}}}
/***
|Name|Plugin: jsMath|
|Created by|BobMcElrath|
|Email|my first name at my last name dot org|
|Location|http://bob.mcelrath.org/tiddlyjsmath.html|
|Version|1.6|
|Requires|[[TiddlyWiki|http://www.tiddlywiki.com]] ≥ 2.0.3, [[jsMath|http://www.math.union.edu/~dpvc/jsMath/]] ≥ 3.0|
!Description
LaTeX is the world standard for specifying, typesetting, and communicating mathematics among scientists, engineers, and mathematicians. For more information about LaTeX itself, visit the [[LaTeX Project|http://www.latex-project.org/]]. This plugin typesets math using [[jsMath|http://www.math.union.edu/~dpvc/jsMath/]], which is an implementation of the TeX math rules and typesetting in javascript, for your browser. Notice the small button in the lower right corner which opens its control panel.
!Installation
In addition to this plugin, you must also [[install jsMath|http://www.math.union.edu/~dpvc/jsMath/download/jsMath.html]] on the same server as your TiddlyWiki html file. If you're using TiddlyWiki without a web server, then the jsMath directory must be placed in the same location as the TiddlyWiki html file.
I also recommend modifying your StyleSheet use serif fonts that are slightly larger than normal, so that the math matches surrounding text, and \\small fonts are not unreadable (as in exponents and subscripts).
{{{
.viewer {
line-height: 125%;
font-family: serif;
font-size: 12pt;
}
}}}
You may also optionally add the following code to load jsMath in [[MarkupPostHead]] if you desire.
{{{
<!--{{{-->
<script src="jsMath/jsMath.js"></script>
<!--}}}-->
}}}
[[Plugin: jsMath]] will normally load jsMath dynamically using AJAX, but adding the above in [[MarkupPostHead]] may be useful if you have jsMath stored in a non-standard location, or if your browser's cross-site origin policy forbids loading files from file URL's using AJAX. (e.g. Google Chrome)
!History
* 11-Nov-05, version 1.0, Initial release
* 22-Jan-06, version 1.1, updated for ~TW2.0, tested with jsMath 3.1, editing tiddlywiki.html by hand is no longer necessary.
* 24-Jan-06, version 1.2, fixes for Safari, Konqueror
* 27-Jan-06, version 1.3, improved error handling, detect if ajax was already defined (used by ZiddlyWiki)
* 12-Jul-06, version 1.4, fixed problem with not finding image fonts
* 26-Feb-07, version 1.5, fixed problem with Mozilla "unterminated character class".
* 27-Feb-07, version 1.5.1, Runs compatibly with TW 2.1.0+, by Bram Chen
* 5-May-11, version 1.6, Use a script tag to load in Chrome, use jQuery for ajax
!Examples
|!Source|!Output|h
|{{{The variable $x$ is real.}}}|The variable $x$ is real.|
|{{{The variable \(y\) is complex.}}}|The variable \(y\) is complex.|
|{{{This \[\int_a^b x = \frac{1}{2}(b^2-a^2)\] is an easy integral.}}}|This \[\int_a^b x = \frac{1}{2}(b^2-a^2)\] is an easy integral.|
|{{{This $$\int_a^b \sin x = -(\cos b - \cos a)$$ is another easy integral.}}}|This $$\int_a^b \sin x = -(\cos b - \cos a)$$ is another easy integral.|
|{{{Block formatted equations may also use the 'equation' environment \begin{equation} \int \tan x = -\ln \cos x \end{equation} }}}|Block formatted equations may also use the 'equation' environment \begin{equation} \int \tan x = -\ln \cos x \end{equation}|
|{{{Equation arrays are also supported \begin{eqnarray} a &=& b \\ c &=& d \end{eqnarray} }}}|Equation arrays are also supported \begin{eqnarray} a &=& b \\ c &=& d \end{eqnarray} |
|{{{I spent \$7.38 on lunch.}}}|I spent \$7.38 on lunch.|
|{{{I had to insert a backslash (\\) into my document}}}|I had to insert a backslash (\\) into my document|
!Code
***/
//{{{
// Load jsMath
if(typeof jsMath == 'undefined') {
jsMath = {
Setup: {inited: 1}, // don't run jsMath.Setup.Body() yet
Autoload: {root: new String(document.location).replace(/[^\/]*$/,'jsMath/')} // URL to jsMath directory, change if necessary
};
try {
jQuery.ajax({url: jsMath.Autoload.root+"jsMath.js", dataType: 'script',
async: false, error: function(j, s, e) { throw(e); } });
jsMath.Setup.inited=0; // allow jsMath.Setup.Body() to run again
} catch(e) {
if(navigator.userAgent.toLowerCase().indexOf('chrome') > -1
&& document.location.toString().indexOf('file') == 0) {
var jsMathScript = '<!--{{{-->\n<script src="jsMath/jsMath.js"></script>\n<!--}}}-->';
var MarkupPostHead = store.fetchTiddler("MarkupPostHead");
if(MarkupPostHead && !MarkupPostHead.match(jsMathScript)) {
MarkupPostHead.text += jsMathScript;
} else {
MarkupPostHead = store.createTiddler("MarkupPostHead");
MarkupPostHead.text = jsMathScript;
}
throw("jsMath added to MarkupPostHead: now save and reload (Google Chrome cross origin file:/// URL workaround)");
} else {
alert("jsMath was not found: you must place the 'jsMath' directory in the same place as this file. "
+"The error was:\n"+e.name+": "+e.message);
throw(e); // abort eval
}
}
}
// Define wikifers for latex
config.formatterHelpers.mathFormatHelper = function(w) {
var e = document.createElement(this.element);
e.className = this.className;
var endRegExp = new RegExp(this.terminator, "mg");
endRegExp.lastIndex = w.matchStart+w.matchLength;
var matched = endRegExp.exec(w.source);
if(matched) {
var txt = w.source.substr(w.matchStart+w.matchLength,
matched.index-w.matchStart-w.matchLength);
if(this.keepdelim) {
txt = w.source.substr(w.matchStart, matched.index+matched[0].length-w.matchStart);
}
e.appendChild(document.createTextNode(txt));
w.output.appendChild(e);
w.nextMatch = endRegExp.lastIndex;
}
}
config.formatters.push({
name: "displayMath1",
match: "\\\$\\\$",
terminator: "\\\$\\\$\\n?", // 2.0 compatability
termRegExp: "\\\$\\\$\\n?",
element: "div",
className: "math",
handler: config.formatterHelpers.mathFormatHelper
});
config.formatters.push({
name: "inlineMath1",
match: "\\\$",
terminator: "\\\$", // 2.0 compatability
termRegExp: "\\\$",
element: "span",
className: "math",
handler: config.formatterHelpers.mathFormatHelper
});
var backslashformatters = new Array(0);
backslashformatters.push({
name: "inlineMath2",
match: "\\\\\\\(",
terminator: "\\\\\\\)", // 2.0 compatability
termRegExp: "\\\\\\\)",
element: "span",
className: "math",
handler: config.formatterHelpers.mathFormatHelper
});
backslashformatters.push({
name: "displayMath2",
match: "\\\\\\\[",
terminator: "\\\\\\\]\\n?", // 2.0 compatability
termRegExp: "\\\\\\\]\\n?",
element: "div",
className: "math",
handler: config.formatterHelpers.mathFormatHelper
});
backslashformatters.push({
name: "displayMath3",
match: "\\\\begin\\{equation\\}",
terminator: "\\\\end\\{equation\\}\\n?", // 2.0 compatability
termRegExp: "\\\\end\\{equation\\}\\n?",
element: "div",
className: "math",
handler: config.formatterHelpers.mathFormatHelper
});
// These can be nested. e.g. \begin{equation} \begin{array}{ccc} \begin{array}{ccc} ...
backslashformatters.push({
name: "displayMath4",
match: "\\\\begin\\{eqnarray\\}",
terminator: "\\\\end\\{eqnarray\\}\\n?", // 2.0 compatability
termRegExp: "\\\\end\\{eqnarray\\}\\n?",
element: "div",
className: "math",
keepdelim: true,
handler: config.formatterHelpers.mathFormatHelper
});
// The escape must come between backslash formatters and regular ones.
// So any latex-like \commands must be added to the beginning of
// backslashformatters here.
backslashformatters.push({
name: "escape",
match: "\\\\.",
handler: function(w) {
w.output.appendChild(document.createTextNode(w.source.substr(w.matchStart+1,1)));
w.nextMatch = w.matchStart+2;
}
});
config.formatters=backslashformatters.concat(config.formatters);
window.wikify = function(source,output,highlightRegExp,tiddler)
{
if(source && source != "") {
if(version.major == 2 && version.minor > 0) {
var wikifier = new Wikifier(source,getParser(tiddler),highlightRegExp,tiddler);
wikifier.subWikifyUnterm(output);
} else {
var wikifier = new Wikifier(source,formatter,highlightRegExp,tiddler);
wikifier.subWikify(output,null);
}
jsMath.ProcessBeforeShowing();
}
}
//}}}
Consider the following dataset: [img(50%+,auto)[http://homepages.spa.umn.edu/~vinals/tspot_files/pca_original.jpg]]
The correlation matrix and related eigenvalues are given by, [img(100%+,auto)[http://homepages.spa.umn.edu/~vinals/tspot_files/pca_algebra.jpg]]
A one dimensional projection can be obtained by the PCA method, yielding [img(50%+,auto)[http://homepages.spa.umn.edu/~vinals/tspot_files/pca_1d.jpg]]
Including the eigenvector with the two largest eigenvalues, one would have a two dimensional projection:
[img(50%+,aut[img(40%+,auto)[http://homepages.spa.umn.edu/~vinals/tspot_files/pca_2d.jpg]] [img(50%+,aut[img(60%+,auto)[http://homepages.spa.umn.edu/~vinals/tspot_files/pca_2d_flat.jpg]]
# Develop the high level structure of the algorithm first, and split it in self contained modules.
# Write and test the modules and subroutines/functions individually. Unless the code is very short, do not try to write the entire program at once.
# Use versioning software and repositories for complex code - else keep safety copies. The University operates its own ~GitHub site https://github.umn.edu/ (use your x500 credentials to log in). Keeping versions allows you to go back to previous versions and restore working copies of your project. You can find help at https://help.github.com/articles/set-up-git and https://help.github.com/articles/create-a-repo.
# Stick as much as possible to standard libraries and versions. Many codes need to be run on several different systems and compatibility is a frequent issue.
# Debugging a particular code is the lengthiest part of the development. A code well writen (modular, consistent in notation) makes the task of finding errors simpler. If modules have been tested independently, even better.
# In debugging, if possible, see if the code (or parts) reproduce known solutions.
# Keep in mind that it is possible (although rare) that the computer may make a mistake (typically though internal code optimization).
# Style:
## Indentation is a great way to be able to read, understand and debug the code.
## Comments are useful to document and return to the code at a later time.
## Introduce consisten naming of variables, indices (i,j, ...), loops and arrays.
!!!Data types
|!Type |!FORTRAN descriptor |!C descriptor |!Data size |
|Character | character | char | 1 byte |
|Integer | integer | int | 4 bytes |
|Real (single precision) | real(4) | float | 4 bytes |
|Real (double precision) | real(8) | double | 8 bytes |
Remember:
* 1 byte = 8 bits
* Modern ~CPUs assume 64 bit words by default. Processor word length is in fact higher. An Intel Haswell processor, for example, has registers that are 512 bit wide (i.e., eight double precision real numbers can be operated on in parallel). This leads to effectively 16 double precision operations per CPU cycle. If interested in details, check http://www.anandtech.com/show/6355/intels-haswell-architecture/8
!!! Lecture 12. October 5, 2020. Numerical integration
* [[Class notes|https://github.umn.edu/vinals/PHYS4041-2020/blob/master/notes/[[Class notes|https://github.umn.edu/vinals/PHYS4041-2020/blob/master/notes/numerical_integration.ipynb]].
!!! Lecture 13. October 6, 2020. Extrapolation
* [[Class notes|https://github.umn.edu/vinals/PHYS4041-2020/blob/master/notes/[[Class notes|https://github.umn.edu/vinals/PHYS4041-2020/blob/master/notes/extrapolation.ipynb]].
!!! Lecture 14. October 7, 2020. Introduction to Ordinary Differential Equations I
* Introdcution, forward and backward Euler method.
* Considerations about error propagation. Local error in the Euler method.
* [[Class notes|https://github.umn.edu/vinals/PHYS4041-2020/blob/master/notes/[[Class notes|https://github.umn.edu/vinals/PHYS4041-2020/blob/master/notes/ODE_Initial_Value_Problems.ipynb]].
!!! Lecture 15. October 12, 2020. Introduction to Ordinary Differential Equations II
* Higher order methods. ~Runge-Kutta methods.
* [[Class notes|https://github.umn.edu/vinals/PHYS4041-2020/blob/master/notes/[[Class notes|https://github.umn.edu/vinals/PHYS4041-2020/blob/master/notes/ODE_Initial_Value_Problems.ipynb]].
!!! Lecture 16. October 13, 2020. Introduction to Ordinary Differential Equations III
* Multi-step methods.
* Implicit algorithms and the ~Predictor-Corrector method.
* Systems of differential equations, and higher order equations.
* [[Class notes|https://github.umn.edu/vinals/PHYS4041-2020/blob/master/notes/[[Class notes|https://github.umn.edu/vinals/PHYS4041-2020/blob/master/notes/ODE_Initial_Value_Problems.ipynb]].
!!! Lecture 17. October 14, 2020. Stability of ODE algorithms.
* [[Class notes|https://github.umn.edu/vinals/PHYS4041-2020/blob/master/notes/[[Class notes|https://github.umn.edu/vinals/PHYS4041-2020/blob/master/notes/ODE_stability.ipynb]].
!!! Lecture 18. October 19, 2020. Boundary Value and Eigenvalue Problems
* [[Class notes|https://github.umn.edu/vinals/PHYS4041-2020/blob/master/notes/[[Class notes|https://github.umn.edu/vinals/PHYS4041-2020/blob/master/notes/ODE_boundary _value.ipynb]].
!!! Lecture 19. October 20, 2020. Introduction to Partial Differential Equations
* [[Class notes|https://github.umn.edu/vinals/PHYS4041-2020/blob/master/notes/[[Class notes|https://github.umn.edu/vinals/PHYS4041-2020/blob/master/notes/PDE_Introduction.ipynb]].
!!! Lecture 20. October 21, 2020. Algorithms for Hyperbolic Partial Differential Equations
* [[Class notes|https://github.umn.edu/vinals/PHYS4041-2020/blob/master/notes/[[Class notes|https://github.umn.edu/vinals/PHYS4041-2020/blob/master/notes/PDE_hyperbolic.ipynb]].
!!! Lecture 1. September 8, 2020. Course introduction
Introduction to PHYS/AST 4041. Online resources and class mechanics.
!!! Lecture 2. September 9, 2020. Number representation on a computer
* [[Class notes|https://github.umn.edu/vinals/PHYS4041-2020/blob/master/notes/number_representation.ipynb]].
* The Institute of Electrical and Electronics Engineers (IEEE) develops and maintains a number of standards for various aspects of computing. In particular, you can read the [[2019 IEEE Standard for Floating-Point Arithmetic|http://homepages.spa.umn.edu/~vinals/tspot_files/phys4041/2020/IEEE Standard 754-2019.pdf]] (provided by Richard Schwarz).
* [Landau]: Section 1.5.
!!! Lecture 3. September 14, 2020. Computer Hardware.
* [[Computer Hardware]]
* [Landau]: Chapter 14. More than you really need to know at this stage, but it is an interesting read.
* General considerations about errors. Error propagation. [[Class notes|https://github.umn.edu/vinals/PHYS4041-2020/blob/master/notes/error_propagation.ipynb]].
!!! Lecture 4. September 15, 2020. Review of Calculus (Taylor series)
* [[Class notes|https://github.umn.edu/vinals/PHYS4041-2020/blob/master/notes/calculus_review.ipynb]].
* Any calculus book that you may have used as an undergraduate. There is also an interesting article in [[Wikipedia|https://en.wikipedia.org/wiki/Taylor_series]].
* As discussed in class, here is the [[IBM Accurate Mathematical Library|https://sourceware.org/git/?p=glibc.git;a=blob;f=sysdeps/ieee754/dbl-64/s_sin.c;hb=HEAD#l281]] code to compute $\sin(x)$ and $\cos(x)$ by use of Taylor series.
!!! Lecture 5. September 16, 2020. Basic ideas of numerical analysis
Iteration, local approximation, and numerical instability.
* [[Class notes|https://github.umn.edu/vinals/PHYS4041-2020/blob/master/notes/basic_ideas_numerical_analysis.ipynb]]
!!! Lecture 6. September 21, 2020. Finding the zeros of a given function $f(x) = 0$
* [[Class notes|https://github.umn.edu/vinals/PHYS4041-2020/blob/master/notes/function_zeros.ipynb]]
* [Landau]: Pages 152-157. We will not cover yet the case in which $f$ is a function of many variables. If you are interested, or need it, it is discussed in great detail in [Landau], pages 160-193.
!!! Lecture 7. September 22, 2020. Linear Systems I: Direct Methods
We will follow [Vesely], Chapter 2. The following topics will be covered.
* Matrix algebra. Statement of the problem: find the solution of a system of equations $A \mathbf{x} = \mathbf{b}$, where the coefficient matrix $A$ is known, as well as the right hand side vector $\mathbf{b}$. The problem is stated as finding the solution vector $\mathbf{x}$.
* Direct methods: backward substitution and Gauss elimination.
* Special matrices, and operation count for the solution of the system. [[Tridiagonal system solver|https://github.umn.edu/vinals/PHYS4041-2020/blob/master/notes/tridiagonal_solution.ipynb]].
!!! Lecture 8. September 23, 2020. Linear Systems II: Error Analysis
Error analysis. Ill conditioned problem.
* [[Class notes|https://github.umn.edu/vinals/PHYS4041-2020/blob/master/notes/condition_number.ipynb]]
!!! Lecture 9. September 28, 2020. Linear Systems III: Iterative methods.
* [[Class notes|https://github.umn.edu/vinals/PHYS4041-2020/blob/master/notes/[[Class notes|https://github.umn.edu/vinals/PHYS4041-2020/blob/master/notes/linear_systems_iterative_methods.ipynb]].
!!! Lecture 10. September 29, 2020. Linear Systems IV: Eigenvalue problems.
* [[Class notes|https://github.umn.edu/vinals/PHYS4041-2020/blob/master/notes/[[Class notes|https://github.umn.edu/vinals/PHYS4041-2020/blob/master/notes/linear_systems_eigenvalue_problems.ipynb]].
!!! Lecture 11. September 30, 2020. Numerical differentiation.
* [[Class notes|https://github.umn.edu/vinals/PHYS4041-2020/blob/master/notes/[[Class notes|https://github.umn.edu/vinals/PHYS4041-2020/blob/master/notes/numerical_differentiation.ipynb]].
University of Minnesota. School of Physics and Astronomy
PHYS/AST 4041. Computational Methods in the Physical Sciences <br><br>
@media print {#mainMenu {display: none ! important;}}
@media print {#topMenu {display: none ! important;}}
@media print {#sidebar {display: none ! important;}}
@media print {#messageArea {display: none ! important;}}
@media print {#toolbar {display: none ! important;}}
@media print {.header {display: none ! important;}}
@media print {.tiddler .subtitle {display: none ! important;}}
@media print {.tiddler .toolbar {display; none ! important; }}
@media print {.tiddler .tagging {display; none ! important; }}
@media print {.tiddler .tagged {display; none ! important; }}
.justifyright {
text-align: right;
}
Supplementary reading material for the topics covered in class.
# [[Basic introduction|http://homepages.spa.umn.edu/~vinals/classes/phys4041/phys4041_notes.pdf]] to scientific programming elements which we will be using in class.
# [[Computer Hardware]].
# [[Programming Tips]].
# C [[code|https://sourceware.org/git/?p=glibc.git;a=blob;f=sysdeps/ieee754/dbl-64/s_sin.c;hb=HEAD#l281]] for computing sin(x) or cos(x).
# [[Central Limit Theorem|https://anydice.com/]]. An interactive demonstration of the limit.
# LU Decomposition (John Burkardt, Florida State University, http://people.sc.fsu.edu/~jburkardt/classes/gateway_2014/gateway_2014.html ). [[pdf|http://homepages.spa.umn.edu/~vinals/tspot_files/phys4041/2014/LUDecomposition.pdf]].
# The Conjugate Gradient [[method|http://www.cs.cmu.edu/~quake-papers/painless-conjugate-gradient.pdf]].
# Examples of PDE solutions (Java applets).
## [[Poisson equation|https://phet.colorado.edu/en/simulation/charges-and-fields]].
## [[Vibrating string|https://phet.colorado.edu/sims/html/wave-on-a-string/latest/wave-on-a-string_en.html]].
## [[Ising model Monte Carlo|https://physics.weber.edu/schroeder/software/demos/IsingModel.html]].
## [[Molecular Dynamics of a fluid|https://physics.weber.edu/schroeder/software/demos/MDv0.html]].
## [[One dimensional, time dependent Schrodinger equation|http://www.math.stonybrook.edu/~ccc/WaveSim.html]].
## [[Assorted examples|http://www.falstad.com/mathphysics.html]]. 2D waves and Poisson.
## [[One dimensional diffusion equation|https://ocw.mit.edu/ans7870/18/18.03/s06/tools/HeatEquation.html]].
## [[A collection of links|http://www.math.lamar.edu/faculty/maesumi/PDE1.html]].
## [[More links to HTML5|http://physics.bu.edu/~duffy/HTML5/]].
# [[Ewald Sums|http://homepages.spa.umn.edu/~vinals/tspot_files/phys4041/2014/EwaldSum.pdf]].
# Molecular Dynamics simulations.
## [[An introduction|http://www.youtube.com/watch?v=lLFEqKl3sm4&list=PLY5jumQrwEXXlVEg1bsUR72yc8v3t4dJc]]
## [[Water|http://www.youtube.com/watch?v=x8Atqz5YvzQ]].
## [[Peptide|http://www.youtube.com/watch?v=WScPbvUwDno]]
## [[Ice melting|https://www.youtube.com/watch?v=NQhjAtCKghE]].
## A simple deterministic and time reversal invariant thermostat. Henk van [[Beijeren|http://arxiv.org/abs/1411.2983]].
# [[Principal Component Analysis example]].
# [[K-means]].
# [[Expectation Maximization Clustering]] example.
# [[Example of pseudo-code|http://homepages.spa.umn.edu/~vinals/tspot_files/phys4041/2018/p_code.pdf]].
[[Additional Reading]], not directly related to the topics covered in class.
!Course Information
''Instructor:'' Jorge Viñals
''Office:'' Tate Hall, 130-17.
''Contact Information:'' vinals@umn.edu
''Office Hours: '' MW, 9:00 a.m. - 10:00 a.m.
''Lecture schedule: '' MTW 8:00 a.m. - 8:50 a.m.
''Lab schedule:'' Sessions 1 and 2: Th 8:00 a.m. - 8:50 a.m. Sessions 3 and 4: Th 9:00 a.m. - 9:50 a.m. Keller Hall 1-250 has been reserved.
''TA:'' Peter Stoeckl (stoec034@umn.edu) and Rafid Mahbub (mahbu004@umn.edu).
''TA Office Hours:'' Peter: M 12:00 p.m. - 1:00 p.m. ([[umn.zoom.us/j/94849792545|https://umn.zoom.us/j/94849792545]]). Rafid: F 2:00 p.m. - 3:00 p.m. ([[umn.zoom.us/j/7695447358|https://umn.zoom.us/j/7695447358]]).
!Required Materials
There is no required textbook for this class. Class notes and material posted online should be sufficient. There are, however, a few textbooks that are available online, at the University Library, and for sale at the University Bookstore [[(2018)|https://www.bookstores.umn.edu/course/phys-4041-lec-001-computational-methods-twin-cities-fall-2018]] and [[(2020)|https://www.bookstores.umn.edu/course/phys-4041-lec-001-computational-methods-twin-cities-fall-2020]] that follow and complement the material in this class.
* F.J. Vesely, "Computational Physics - An Introduction". Kluwer Academic/Plenum Press, New York 2001, second edition. ISBN: 0-306-46631-7. This book is available [[online|https://link-springer-com.ezp2.lib.umn.edu/book/10.1007%2F978-1-4757-2307-6]] for University students.
* R.H. Landau, M. Paez, and C. Bordeianu, "A Survey of Computational Physics". Princeton University Press, New Jersey, 2008. ISBN: 978-0-691-13137-5. [[Online access| https://www-jstor-org.ezp3.lib.umn.edu/stable/j.ctv7h0rx7]]. Also at Walter Library (~QC52.L36 2008).
* B.A. Stickler and E. Schachinger. "Basic Concepts in Computational Physics". Springer, New York (2014), second edition. This book is available [[online|http://link.springer.com/book/10.1007%2F978-3-319-02435-6]] for University students.
A generally good condensed reference, with example codes is:
* W.H. Press, S.A. Teukolsky, W.T. Vetterling, and B.P. Flannery, [[Numerical recipes|http://numerical.recipes/]]: The Art of Scientific Computing. Cambridge University Press, (2007) third edition.
!Course goals and objectives
This class is an introduction to Scientific Computing from the point of view of Physics applications. It is not a programming class, although some programming experience is assumed and will be necessary during the class. This class is not about learning to use algorithms written by others either, although some reuse, especially of general purpose libraries, will be essential. Rather, the class will focus on the understanding and design of a computing strategy to solve physics problems, the basic numerical analysis that is necessary to develop classes of algorithms, and the implementation of simple but custom application codes.
!!!Topics
# Introduction to scientific computing. Computer architecture, number representation and approximations. Errors and propagation.
# Basic elements of numerical analysis: iteration, local approximation, convergence and accuracy, and stability.
# Brief survey of classical problems: linear and nonlinear systems, interpolation and extrapolation. Differentiation, integration and solution of systems or ordinary differential equations. Boundary and eigenvalue problems.
# Partial differential equations. Methods for elliptic, parabolic and hyperbolic equations. Relaxation methods for boundary value problems.
# Statistical description of data. Probability distributions and their estimation. Data modeling and regression.
# Monte Carlo methods and their quantum extension.
# Elements of data mining. Classification and decision trees. Neural networks. Association analysis in Big Data. Cluster analysis.
# A high level view of more modern architectures: many core chips, and tightly coupled clusters. Multi threaded and distributed memory programming.
! Homeworks, Laboratory sessions, Exams, and Grading
* __Homeworks__ will be assigned several times during the semester. Typically, the homework will be due one week after being assigned. There will be reduced credit for late homeworks: 75% for up to one week late, 50% for up to two weeks late. No credit after that. Homeworks can be completed either individually, or in a team of two people.
* __A "laboratory session"__ will be held each week. Assume that the Lab session will be conducted remotely. Four different lab sessions will be held each week to allow work in small groups. Each student is expected to attend only one of the weekly sessions.
* __Mid term and final exams__ will the "take home" exams. The structure of the exams is as follows: One long term, programming assignment will be given to each of you. The mid term exam will cover a description of the physical problem, defining mathematical equations (if appropriate), algorithm, and any accuracy or stability considerations. ''A written report covering these aspects of the problem is due on October 31, 2020''. The final exam will require the completed code, and numerical results from its execution. Depending on the individual case, an oral presentation may be required. ''A final report containing the description of the algorithm, the file(s) containing the code, and numerical results as appropriate are due on December 17, 2020''.
* __Grading__. Homeworks: 40%. Mid Term Exam: 10%, Final Exam: 35%, Laboratory and general class engagement: 15%.
!Attendance policy
Class attendance is required.
!Departamental policies
''ATHLETES'' must provide their official University of Minnesota athletic letter containing the approved competition schedule to their instructor and the staff in office 148. Away exams will be arranged with the athletic adviser traveling with the team. Accommodations will be made for official university sports only (i.e. no accommodations will be made for intramurals, club sports, etc.)
''DISABILITY SERVICES:'' If you have accommodations for this course, please provide the staff in the office (Tate 130) with a copy of your accommodation letter for the current semester. Exams will be arranged according to accommodations and sent to the testing center for administration.
!Mandatory policy information
!!!Student Conduct Code
The University seeks an environment that promotes academic achievement and integrity, that is protective of free inquiry, and that serves the educational mission of the University. Similarly, the University seeks a community that is free from violence, threats, and intimidation; that is respectful of the rights, opportunities, and welfare of students, faculty, staff, and guests of the University; and that does not threaten the physical or mental health or safety of members of the University community. As a student at the University you are expected adhere to Board of Regents Policy: Student Conduct Code. To review the Student Conduct Code, please see: [[Regents Student Conduct Code|http://regents.umn.edu/sites/default/files/policies/Student_Conduct_Code.pdf]]
Note that the conduct code specifically addresses disruptive classroom conduct, which means "engaging in behavior that substantially or repeatedly interrupts either the instructor's ability to teach or student learning. The classroom extends to any setting where a student is engaged in work toward academic credit or satisfaction of program-based requirements or related activities."
!!!Scholastic dishonesty
You are expected to do your own academic work and cite sources as necessary. Failing to do so is scholastic dishonesty. Scholastic dishonesty means plagiarizing; cheating on assignments or examinations; engaging in unauthorized collaboration on academic work; taking, acquiring, or using test materials without faculty permission; submitting false or incomplete records of academic achievement; acting alone or in cooperation with another to falsify records or to obtain dishonestly grades, honors, awards, or professional endorsement; altering, forging, or misusing a University academic record; or fabricating or falsifying data, research procedures, or data analysis ([[Student Conduct Code|http://regents.umn.edu/sites/default/files/policies/Student_Conduct_Code.pdf]]).
If it is determined that a student has cheated, he or she may be given an "F" or an "N" for the course, and may face additional sanctions from the University. For additional information, please see: http://policy.umn.edu/Policies/Education/Education/INSTRUCTORRESP.html.
The Office for Student Conduct and Academic Integrity has compiled a useful list of Frequently Asked Questions pertaining to scholastic dishonesty: http://www1.umn.edu/oscai/integrity/student/index.html. If you have additional questions, please clarify with your instructor for the course. Your instructor can respond to your specific questions regarding what would constitute scholastic dishonesty in the context of a particular class-e.g., whether collaboration on assignments is permitted, requirements and methods for citing sources, if electronic aids are permitted or prohibited during an exam.
!!!Disability Accommodations
The University is committed to providing quality education to all students regardless of ability. Determining appropriate disability accommodations is a collaborative process. You as a student must register with Disability Services and provide documentation of your disability. The course instructor must provide information regarding a course's content, methods, and essential components. The combination of this information will be used by Disability Services to determine appropriate accommodations for a particular student in a particular course. For more information, please reference Disability Services: https://diversity.umn.edu/disability/.
If you have a letter from the Disability Resource Center regarding a special accommodation, please make sure that you share it with the instructor, and the Physics front office.
!!!Use of Personal Electronic Devices in the Classroom
Using personal electronic devices in the classroom setting can hinder instruction and learning, not only for the student using the device but also for other students in the class. To this end, the University establishes the right of each faculty member to determine if and how personal electronic devices are allowed to be used in the classroom. For complete information, please reference: http://www.policy.umn.edu/Policies/Education/Education/STUDENTRESP.html
!!!Makeup Work for Legitimate Absences
Students will not be penalized for absence during the semester due to unavoidable or legitimate circumstances. Such circumstances include verified illness, participation in intercollegiate athletic events, subpoenas, jury duty, military service, bereavement, and religious observances. Such circumstances do not include voting in local, state, or national elections. For complete information, please see: http://policy.umn.edu/Policies/Education/Education/MAKEUPWORK.html.
!!!Appropriate Student Use of Class Notes and Course Materials
Taking notes is a means of recording information but more importantly of personally absorbing and integrating the educational experience. However, broadly disseminating class notes beyond the classroom community or accepting compensation for taking and distributing classroom notes undermines instructor interests in their intellectual work product while not substantially furthering instructor and student interests in effective learning. Such actions violate shared norms and standards of the academic community. For additional information, please see: http://policy.umn.edu/education/studentresp .
!!!Grading and Transcripts
The University utilizes plus and minus grading on a 4.000 cumulative grade point scale in accordance with the following:
|A |4.000 - Represents achievement that is outstanding relative to the level necessary to meet course requirements|
|A- |3.667 |
|B+ |3.333 |
|B |3.000 - Represents achievement that is significantly above the level necessary to meet course requirements |
|B- |2.667 |
|C+ |2.333 |
|C |2.000 - Represents achievement that meets the course requirements in every respect |
|C- |1.667 |
|D+ |1.333 |
|D |1.000 - Represents achievement that is worthy of credit even though it fails to meet fully the course requirements |
|S |Represents achievement that is satisfactory, which is equivalent to a C- or better. |
For additional information, please refer to: http://policy.umn.edu/Policies/Education/Education/GRADINGTRANSCRIPTS.html.
!!!Sexual Harassment
"Sexual harassment" means unwelcome sexual advances, requests for sexual favors, and/or other verbal or physical conduct of a sexual nature. Such conduct has the purpose or effect of unreasonably interfering with an individual's work or academic performance or creating an intimidating, hostile, or offensive working or academic environment in any University activity or program. Such behavior is not acceptable in the University setting. For additional information, please consult Board of Regents Policy: https://policy.umn.edu/hr/sexharassassault .
!!!Equity, Diversity, Equal Opportunity, and Affirmative Action
The University will provide equal access to and opportunity in its programs and facilities, without regard to race, color, creed, religion, national origin, gender, age, marital status, disability, public assistance status, veteran status, sexual orientation, gender identity, or gender expression. For more information, please consult Board of Regents Policy: http://regents.umn.edu/sites/regents.umn.edu/files/policies/Equity_Diversity_EO_AA.pdf .
!!!Mental Health and Stress Management
As a student you may experience a range of issues that can cause barriers to learning, such as strained relationships, increased anxiety, alcohol/drug problems, feeling down, difficulty concentrating and/or lack of motivation. These mental health concerns or stressful events may lead to diminished academic performance and may reduce your ability to participate in daily activities. University of Minnesota services are available to assist you. You can learn more about the broad range of confidential mental health services available on campus via the Student Mental Health Website: http://www.mentalhealth.umn.edu.
!!! Academic Freedom and Responsibility: for courses that do not involve students in research
Academic freedom is a cornerstone of the University. Within the scope and content of the course as defined by the instructor, it includes the freedom to discuss relevant matters in the classroom. Along with this freedom comes responsibility. Students are encouraged to develop the capacity for critical judgment and to engage in a sustained and independent search for truth. Students are free to take reasoned exception to the views offered in any course of study and to reserve judgment about matters of opinion, but they are responsible for learning the content of any course of study for which they are enrolled.
Reports of concerns about academic freedom are taken seriously, and there are individuals and offices available for help. Contact the instructor, the Department Chair, your adviser, the associate dean of the college, or the Vice Provost for Faculty and Academic Affairs in the Office of the Provost.
<<cloud systemConfig systemTiddlers DiscoveryPackage IconPackage NavigationPackage>>
/***
|Name|TagCloudPlugin|
|Source|http://www.TiddlyTools.com/#TagCloudPlugin|
|Version|1.7.0|
|Author|Eric Shulman|
|Original Author|Clint Checketts|
|License|http://www.TiddlyTools.com/#LegalStatements|
|~CoreVersion|2.1|
|Type|plugin|
|Description|present a 'cloud' of tags (or links) using proportional font display|
!Usage
<<<
{{{
<<cloud type action:... limit:... tag tag tag ...>>
<<cloud type action:... limit:... +TiddlerName>>
<<cloud type action:... limit:... -TiddlerName>>
<<cloud type action:... limit:... =tagvalue>>
}}}
where:
* //type// is a keyword, one of:
** ''tags'' (default) - displays a cloud of tags, based on frequency of use
** ''links'' - displays a cloud of tiddlers, based on number of links //from// each tiddler
** ''references'' - displays a cloud of tiddlers, based on number of links //to// each tiddler
* ''action:popup'' (default) - clicking a cloud item shows a popup with links to related tiddlers<br>//or//<br> ''action:goto'' - clicking a cloud item immediately opens the tiddler corresponding to that item
* ''limit:N'' (optional) - restricts the cloud display to only show the N most popular tags/links
* ''tag tag tag...'' (or ''title title title'' if ''links''/''references'' is used)<br>shows all tags/links in the document //except// for those listed as macro parameters
* ''+TiddlerName''<br>show only tags/links read from a space-separated, bracketed list stored in a separate tiddler.
* ''-TiddlerName''<br>show all tags/links //except// those read from a space-separated, bracketed list stored in a separate tiddler.
* ''=tagvalue'' (//only if type=''tags''//)<br>shows only tags that are themselves tagged with the indicated tag value (i.e., ~TagglyTagging usage)
//note: for backward-compatibility, you can also use the macro {{{<<tagCloud ...>>}}} in place of {{{<<cloud ...>>}}}//
<<<
!Examples
<<<
//all tags excluding<<tag systemConfig>>, <<tag excludeMissing>> and <<tag script>>//
{{{<<cloud systemConfig excludeMissing script>>}}}
{{groupbox{<<cloud systemConfig excludeMissing script>>}}}
//top 10 tags excluding<<tag systemConfig>>, <<tag excludeMissing>> and <<tag script>>//
{{{<<cloud limit:10 systemConfig excludeMissing script>>}}}
{{groupbox{<<cloud limit:10 systemConfig excludeMissing script>>}}}
//tags listed in// [[FavoriteTags]]
{{{<<cloud +FavoriteTags>>}}}
{{groupbox{<<cloud +FavoriteTags>>}}}
//tags NOT listed in// [[FavoriteTags]]
{{{<<cloud -FavoriteTags>>}}}
{{groupbox{<<cloud -FavoriteTags>>}}}
//links to tiddlers tagged with 'package'//
{{{<<cloud action:goto =package>>}}}
{{groupbox{<<cloud action:goto =package>>}}}
//top 20 most referenced tiddlers//
{{{<<cloud references limit:20>>}}}
{{groupbox{<<cloud references limit:20>>}}}
//top 20 tiddlers that contain the most links//
{{{<<cloud links limit:20>>}}}
{{groupbox{<<cloud links limit:20>>}}}
<<<
!Revisions
<<<
2009.07.17 [1.7.0] added {{{-TiddlerName}}} parameter to exclude tags that are listed in the indicated tiddler
2009.02.26 [1.6.0] added {{{action:...}}} parameter to apply popup vs. goto action when clicking cloud items
2009.02.05 [1.5.0] added ability to show links or back-links (references) instead of tags and renamed macro to {{{<<cloud>>}}} to reflect more generalized usage.
2008.12.16 [1.4.2] corrected group calculation to prevent 'group=0' error
2008.12.16 [1.4.1] revised tag filtering so excluded tags don't affect calculations
2008.12.15 [1.4.0] added {{{limit:...}}} parameter to restrict the number of tags displayed to the top N most popular
2008.11.15 [1.3.0] added {{{+TiddlerName}}} parameter to include only tags that are listed in the indicated tiddler
2008.09.05 [1.2.0] added '=tagname' parameter to include only tags that are themselves tagged with the specified value (i.e., ~TagglyTagging usage)
2008.07.03 [1.1.0] added 'segments' property to macro object. Extensive code cleanup
<<<
!Code
***/
//{{{
version.extensions.TagCloudPlugin= {major: 1, minor: 7 , revision: 0, date: new Date(2009,7,17)};
//Originally created by Clint Checketts, contributions by Jonny Leroy and Eric Shulman
//Currently maintained and enhanced by Eric Shulman
//}}}
//{{{
config.macros.cloud = {
tagstip: "%1 tiddlers tagged with '%0'",
refslabel: " (%0 references)",
refstip: "%1 tiddlers have links to '%0'",
linkslabel: " (%0 links)",
linkstip: "'%0' has links to %1 other tiddlers",
groups: 9,
init: function() {
config.macros.tagCloud=config.macros.cloud; // for backward-compatibility
config.shadowTiddlers.TagCloud='<<cloud>>';
config.shadowTiddlers.StyleSheetTagCloud=
'/*{{{*/\n'
+'.tagCloud span {line-height: 3.5em; margin:3px;}\n'
+'.tagCloud1{font-size: 80%;}\n'
+'.tagCloud2{font-size: 100%;}\n'
+'.tagCloud3{font-size: 120%;}\n'
+'.tagCloud4{font-size: 140%;}\n'
+'.tagCloud5{font-size: 160%;}\n'
+'.tagCloud6{font-size: 180%;}\n'
+'.tagCloud7{font-size: 200%;}\n'
+'.tagCloud8{font-size: 220%;}\n'
+'.tagCloud9{font-size: 240%;}\n'
+'/*}}}*/\n';
setStylesheet(store.getTiddlerText('StyleSheetTagCloud'),'tagCloudsStyles');
},
getLinks: function(tiddler) { // get list of links to existing tiddlers and shadows
if (!tiddler.linksUpdated) tiddler.changed();
var list=[]; for (var i=0; i<tiddler.links.length; i++) {
var title=tiddler.links[i];
if (store.isShadowTiddler(title)||store.tiddlerExists(title))
list.push(title);
}
return list;
},
handler: function(place,macroName,params) {
// unpack params
var inc=[]; var ex=[]; var limit=0; var action='popup';
var links=(params[0]&¶ms[0].toLowerCase()=='links'); if (links) params.shift();
var refs=(params[0]&¶ms[0].toLowerCase()=='references'); if (refs) params.shift();
if (params[0]&¶ms[0].substr(0,7).toLowerCase()=='action:')
action=params.shift().substr(7).toLowerCase();
if (params[0]&¶ms[0].substr(0,6).toLowerCase()=='limit:')
limit=parseInt(params.shift().substr(6));
while (params.length) {
if (params[0].substr(0,1)=='+') { // read taglist from tiddler
inc=inc.concat(store.getTiddlerText(params[0].substr(1),'').readBracketedList());
} else if (params[0].substr(0,1)=='-') { // exclude taglist from tiddler
ex=ex.concat(store.getTiddlerText(params[0].substr(1),'').readBracketedList());
} else if (params[0].substr(0,1)=='=') { // get tag list using tagged tags
var tagged=store.getTaggedTiddlers(params[0].substr(1));
for (var t=0; t<tagged.length; t++) inc.push(tagged[t].title);
} else ex.push(params[0]); // exclude params
params.shift();
}
// get all items, include/exclude specific items
var items=[];
var list=(links||refs)?store.getTiddlers('title','excludeLists'):store.getTags();
for (var t=0; t<list.length; t++) {
var title=(links||refs)?list[t].title:list[t][0];
if (links) var count=this.getLinks(list[t]).length;
else if (refs) var count=store.getReferringTiddlers(title).length;
else var count=list[t][1];
if ((!inc.length||inc.contains(title))&&(!ex.length||!ex.contains(title)))
items.push({ title:title, count:count });
}
if(!items.length) return;
// sort by decending count, limit results (optional)
items=items.sort(function(a,b){return(a.count==b.count)?0:(a.count>b.count?-1:1);});
while (limit && items.length>limit) items.pop();
// find min/max and group size
var most=items[0].count;
var least=items[items.length-1].count;
var groupSize=(most-least+1)/this.groups;
// sort by title and draw the cloud of items
items=items.sort(function(a,b){return(a.title==b.title)?0:(a.title>b.title?1:-1);});
var cloudWrapper = createTiddlyElement(place,'div',null,'tagCloud',null);
for (var t=0; t<items.length; t++) {
cloudWrapper.appendChild(document.createTextNode(' '));
var group=Math.ceil((items[t].count-least)/groupSize)||1;
var className='tagCloudtag tagCloud'+group;
var tip=refs?this.refstip:links?this.linkstip:this.tagstip;
tip=tip.format([items[t].title,items[t].count]);
if (action=='goto') { // TAG/LINK/REFERENCES GOTO
var btn=createTiddlyLink(cloudWrapper,items[t].title,true,className);
btn.title=tip;
btn.style.fontWeight='normal';
} else if (!links&&!refs) { // TAG POPUP
var btn=createTiddlyButton(cloudWrapper,items[t].title,tip,onClickTag,className);
btn.setAttribute('tag',items[t].title);
} else { // LINK/REFERENCES POPUP
var btn=createTiddlyButton(cloudWrapper,items[t].title,tip,
function(ev) { var e=ev||window.event; var cmt=config.macros.cloud;
var popup = Popup.create(this);
var title = this.getAttribute('tiddler');
var count = this.getAttribute('count');
var refs = this.getAttribute('refs')=='T';
var links = this.getAttribute('links')=='T';
var label = (refs?cmt.refslabel:cmt.linkslabel).format([count]);
createTiddlyLink(popup,title,true);
createTiddlyText(popup,label);
createTiddlyElement(popup,'hr');
if (refs) {
popup.setAttribute('tiddler',title);
config.commands.references.handlePopup(popup,title);
}
if (links) {
var tiddler = store.fetchTiddler(title);
var links=config.macros.cloud.getLinks(tiddler);
for(var i=0;i<links.length;i++)
createTiddlyLink(createTiddlyElement(popup,'li'),
links[i],true);
}
Popup.show();
e.cancelBubble=true; if(e.stopPropagation) e.stopPropagation();
return false;
}, className);
btn.setAttribute('tiddler',items[t].title);
btn.setAttribute('count',items[t].count);
btn.setAttribute('refs',refs?'T':'F');
btn.setAttribute('links',links?'T':'F');
btn.title=tip;
}
}
}
};
//}}}
/***
Description: Contains the stuff you need to use Tiddlyspot
Note, you also need UploadPlugin, PasswordOptionPlugin and LoadRemoteFileThroughProxy
from http://tiddlywiki.bidix.info for a complete working Tiddlyspot site.
***/
//{{{
// edit this if you are migrating sites or retrofitting an existing TW
config.tiddlyspotSiteId = 'phys4041';
// make it so you can by default see edit controls via http
config.options.chkHttpReadOnly = false;
window.readOnly = false; // make sure of it (for tw 2.2)
window.showBackstage = true; // show backstage too
// disable autosave in d3
if (window.location.protocol != "file:")
config.options.chkGTDLazyAutoSave = false;
// tweak shadow tiddlers to add upload button, password entry box etc
with (config.shadowTiddlers) {
SiteUrl = 'http://'+config.tiddlyspotSiteId+'.tiddlyspot.com';
SideBarOptions = SideBarOptions.replace(/(<<saveChanges>>)/,"$1<<tiddler TspotSidebar>>");
OptionsPanel = OptionsPanel.replace(/^/,"<<tiddler TspotOptions>>");
DefaultTiddlers = DefaultTiddlers.replace(/^/,"[[WelcomeToTiddlyspot]] ");
MainMenu = MainMenu.replace(/^/,"[[WelcomeToTiddlyspot]] ");
}
// create some shadow tiddler content
merge(config.shadowTiddlers,{
'TspotControls':[
"| tiddlyspot password:|<<option pasUploadPassword>>|",
"| site management:|<<upload http://" + config.tiddlyspotSiteId + ".tiddlyspot.com/store.cgi index.html . . " + config.tiddlyspotSiteId + ">>//(requires tiddlyspot password)//<br>[[control panel|http://" + config.tiddlyspotSiteId + ".tiddlyspot.com/controlpanel]], [[download (go offline)|http://" + config.tiddlyspotSiteId + ".tiddlyspot.com/download]]|",
"| links:|[[tiddlyspot.com|http://tiddlyspot.com/]], [[FAQs|http://faq.tiddlyspot.com/]], [[blog|http://tiddlyspot.blogspot.com/]], email [[support|mailto:support@tiddlyspot.com]] & [[feedback|mailto:feedback@tiddlyspot.com]], [[donate|http://tiddlyspot.com/?page=donate]]|"
].join("\n"),
'TspotOptions':[
"tiddlyspot password:",
"<<option pasUploadPassword>>",
""
].join("\n"),
'TspotSidebar':[
"<<upload http://" + config.tiddlyspotSiteId + ".tiddlyspot.com/store.cgi index.html . . " + config.tiddlyspotSiteId + ">><html><a href='http://" + config.tiddlyspotSiteId + ".tiddlyspot.com/download' class='button'>download</a></html>"
].join("\n"),
'WelcomeToTiddlyspot':[
"This document is a ~TiddlyWiki from tiddlyspot.com. A ~TiddlyWiki is an electronic notebook that is great for managing todo lists, personal information, and all sorts of things.",
"",
"@@font-weight:bold;font-size:1.3em;color:#444; //What now?// @@ Before you can save any changes, you need to enter your password in the form below. Then configure privacy and other site settings at your [[control panel|http://" + config.tiddlyspotSiteId + ".tiddlyspot.com/controlpanel]] (your control panel username is //" + config.tiddlyspotSiteId + "//).",
"<<tiddler TspotControls>>",
"See also GettingStarted.",
"",
"@@font-weight:bold;font-size:1.3em;color:#444; //Working online// @@ You can edit this ~TiddlyWiki right now, and save your changes using the \"save to web\" button in the column on the right.",
"",
"@@font-weight:bold;font-size:1.3em;color:#444; //Working offline// @@ A fully functioning copy of this ~TiddlyWiki can be saved onto your hard drive or USB stick. You can make changes and save them locally without being connected to the Internet. When you're ready to sync up again, just click \"upload\" and your ~TiddlyWiki will be saved back to tiddlyspot.com.",
"",
"@@font-weight:bold;font-size:1.3em;color:#444; //Help!// @@ Find out more about ~TiddlyWiki at [[TiddlyWiki.com|http://tiddlywiki.com]]. Also visit [[TiddlyWiki.org|http://tiddlywiki.org]] for documentation on learning and using ~TiddlyWiki. New users are especially welcome on the [[TiddlyWiki mailing list|http://groups.google.com/group/TiddlyWiki]], which is an excellent place to ask questions and get help. If you have a tiddlyspot related problem email [[tiddlyspot support|mailto:support@tiddlyspot.com]].",
"",
"@@font-weight:bold;font-size:1.3em;color:#444; //Enjoy :)// @@ We hope you like using your tiddlyspot.com site. Please email [[feedback@tiddlyspot.com|mailto:feedback@tiddlyspot.com]] with any comments or suggestions."
].join("\n")
});
//}}}
| !date | !user | !location | !storeUrl | !uploadDir | !toFilename | !backupdir | !origin |
| 13/10/2020 09:25:11 | Jorge | [[/|http://phys4041.tiddlyspot.com/]] | [[store.cgi|http://phys4041.tiddlyspot.com/store.cgi]] | . | [[index.html | http://phys4041.tiddlyspot.com/index.html]] | . |
| 13/10/2020 13:53:59 | Jorge | [[/|http://phys4041.tiddlyspot.com/]] | [[store.cgi|http://phys4041.tiddlyspot.com/store.cgi]] | . | [[index.html | http://phys4041.tiddlyspot.com/index.html]] | . | ok |
| 13/10/2020 15:02:17 | Jorge | [[/|http://phys4041.tiddlyspot.com/]] | [[store.cgi|http://phys4041.tiddlyspot.com/store.cgi]] | . | [[index.html | http://phys4041.tiddlyspot.com/index.html]] | . | ok |
| 13/10/2020 15:35:35 | Jorge | [[/|http://phys4041.tiddlyspot.com/]] | [[store.cgi|http://phys4041.tiddlyspot.com/store.cgi]] | . | [[index.html | http://phys4041.tiddlyspot.com/index.html]] | . |
| 15/10/2020 08:05:23 | Jorge | [[/|http://phys4041.tiddlyspot.com/]] | [[store.cgi|http://phys4041.tiddlyspot.com/store.cgi]] | . | [[index.html | http://phys4041.tiddlyspot.com/index.html]] | . | failed |
| 15/10/2020 08:05:38 | Jorge | [[/|http://phys4041.tiddlyspot.com/]] | [[store.cgi|http://phys4041.tiddlyspot.com/store.cgi]] | . | [[index.html | http://phys4041.tiddlyspot.com/index.html]] | . | ok |
| 15/10/2020 09:21:34 | Jorge | [[/|http://phys4041.tiddlyspot.com/]] | [[store.cgi|http://phys4041.tiddlyspot.com/store.cgi]] | . | [[index.html | http://phys4041.tiddlyspot.com/index.html]] | . |
| 20/10/2020 14:58:06 | Jorge | [[/|http://phys4041.tiddlyspot.com/]] | [[store.cgi|http://phys4041.tiddlyspot.com/store.cgi]] | . | [[index.html | http://phys4041.tiddlyspot.com/index.html]] | . | ok |
| 20/10/2020 16:02:36 | Jorge | [[/|http://phys4041.tiddlyspot.com/]] | [[store.cgi|http://phys4041.tiddlyspot.com/store.cgi]] | . | [[index.html | http://phys4041.tiddlyspot.com/index.html]] | . | ok |
| 20/10/2020 16:07:58 | Jorge | [[/|http://phys4041.tiddlyspot.com/]] | [[store.cgi|http://phys4041.tiddlyspot.com/store.cgi]] | . | [[index.html | http://phys4041.tiddlyspot.com/index.html]] | . |
/***
|''Name:''|UploadPlugin|
|''Description:''|Save to web a TiddlyWiki|
|''Version:''|4.1.3|
|''Date:''|Feb 24, 2008|
|''Source:''|http://tiddlywiki.bidix.info/#UploadPlugin|
|''Documentation:''|http://tiddlywiki.bidix.info/#UploadPluginDoc|
|''Author:''|BidiX (BidiX (at) bidix (dot) info)|
|''License:''|[[BSD open source license|http://tiddlywiki.bidix.info/#%5B%5BBSD%20open%20source%20license%5D%5D ]]|
|''~CoreVersion:''|2.2.0|
|''Requires:''|PasswordOptionPlugin|
***/
//{{{
version.extensions.UploadPlugin = {
major: 4, minor: 1, revision: 3,
date: new Date("Feb 24, 2008"),
source: 'http://tiddlywiki.bidix.info/#UploadPlugin',
author: 'BidiX (BidiX (at) bidix (dot) info',
coreVersion: '2.2.0'
};
//
// Environment
//
if (!window.bidix) window.bidix = {}; // bidix namespace
bidix.debugMode = false; // true to activate both in Plugin and UploadService
//
// Upload Macro
//
config.macros.upload = {
// default values
defaultBackupDir: '', //no backup
defaultStoreScript: "store.php",
defaultToFilename: "index.html",
defaultUploadDir: ".",
authenticateUser: true // UploadService Authenticate User
};
config.macros.upload.label = {
promptOption: "Save and Upload this TiddlyWiki with UploadOptions",
promptParamMacro: "Save and Upload this TiddlyWiki in %0",
saveLabel: "save to web",
saveToDisk: "save to disk",
uploadLabel: "upload"
};
config.macros.upload.messages = {
noStoreUrl: "No store URL in parmeters or options",
usernameOrPasswordMissing: "Username or password missing"
};
config.macros.upload.handler = function(place,macroName,params) {
if (readOnly)
return;
var label;
if (document.location.toString().substr(0,4) == "http")
label = this.label.saveLabel;
else
label = this.label.uploadLabel;
var prompt;
if (params[0]) {
prompt = this.label.promptParamMacro.toString().format([this.destFile(params[0],
(params[1] ? params[1]:bidix.basename(window.location.toString())), params[3])]);
} else {
prompt = this.label.promptOption;
}
createTiddlyButton(place, label, prompt, function() {config.macros.upload.action(params);}, null, null, this.accessKey);
};
config.macros.upload.action = function(params)
{
// for missing macro parameter set value from options
if (!params) params = {};
var storeUrl = params[0] ? params[0] : config.options.txtUploadStoreUrl;
var toFilename = params[1] ? params[1] : config.options.txtUploadFilename;
var backupDir = params[2] ? params[2] : config.options.txtUploadBackupDir;
var uploadDir = params[3] ? params[3] : config.options.txtUploadDir;
var username = params[4] ? params[4] : config.options.txtUploadUserName;
var password = config.options.pasUploadPassword; // for security reason no password as macro parameter
// for still missing parameter set default value
if ((!storeUrl) && (document.location.toString().substr(0,4) == "http"))
storeUrl = bidix.dirname(document.location.toString())+'/'+config.macros.upload.defaultStoreScript;
if (storeUrl.substr(0,4) != "http")
storeUrl = bidix.dirname(document.location.toString()) +'/'+ storeUrl;
if (!toFilename)
toFilename = bidix.basename(window.location.toString());
if (!toFilename)
toFilename = config.macros.upload.defaultToFilename;
if (!uploadDir)
uploadDir = config.macros.upload.defaultUploadDir;
if (!backupDir)
backupDir = config.macros.upload.defaultBackupDir;
// report error if still missing
if (!storeUrl) {
alert(config.macros.upload.messages.noStoreUrl);
clearMessage();
return false;
}
if (config.macros.upload.authenticateUser && (!username || !password)) {
alert(config.macros.upload.messages.usernameOrPasswordMissing);
clearMessage();
return false;
}
bidix.upload.uploadChanges(false,null,storeUrl, toFilename, uploadDir, backupDir, username, password);
return false;
};
config.macros.upload.destFile = function(storeUrl, toFilename, uploadDir)
{
if (!storeUrl)
return null;
var dest = bidix.dirname(storeUrl);
if (uploadDir && uploadDir != '.')
dest = dest + '/' + uploadDir;
dest = dest + '/' + toFilename;
return dest;
};
//
// uploadOptions Macro
//
config.macros.uploadOptions = {
handler: function(place,macroName,params) {
var wizard = new Wizard();
wizard.createWizard(place,this.wizardTitle);
wizard.addStep(this.step1Title,this.step1Html);
var markList = wizard.getElement("markList");
var listWrapper = document.createElement("div");
markList.parentNode.insertBefore(listWrapper,markList);
wizard.setValue("listWrapper",listWrapper);
this.refreshOptions(listWrapper,false);
var uploadCaption;
if (document.location.toString().substr(0,4) == "http")
uploadCaption = config.macros.upload.label.saveLabel;
else
uploadCaption = config.macros.upload.label.uploadLabel;
wizard.setButtons([
{caption: uploadCaption, tooltip: config.macros.upload.label.promptOption,
onClick: config.macros.upload.action},
{caption: this.cancelButton, tooltip: this.cancelButtonPrompt, onClick: this.onCancel}
]);
},
options: [
"txtUploadUserName",
"pasUploadPassword",
"txtUploadStoreUrl",
"txtUploadDir",
"txtUploadFilename",
"txtUploadBackupDir",
"chkUploadLog",
"txtUploadLogMaxLine"
],
refreshOptions: function(listWrapper) {
var opts = [];
for(i=0; i<this.options.length; i++) {
var opt = {};
opts.push();
opt.option = "";
n = this.options[i];
opt.name = n;
opt.lowlight = !config.optionsDesc[n];
opt.description = opt.lowlight ? this.unknownDescription : config.optionsDesc[n];
opts.push(opt);
}
var listview = ListView.create(listWrapper,opts,this.listViewTemplate);
for(n=0; n<opts.length; n++) {
var type = opts[n].name.substr(0,3);
var h = config.macros.option.types[type];
if (h && h.create) {
h.create(opts[n].colElements['option'],type,opts[n].name,opts[n].name,"no");
}
}
},
onCancel: function(e)
{
backstage.switchTab(null);
return false;
},
wizardTitle: "Upload with options",
step1Title: "These options are saved in cookies in your browser",
step1Html: "<input type='hidden' name='markList'></input><br>",
cancelButton: "Cancel",
cancelButtonPrompt: "Cancel prompt",
listViewTemplate: {
columns: [
{name: 'Description', field: 'description', title: "Description", type: 'WikiText'},
{name: 'Option', field: 'option', title: "Option", type: 'String'},
{name: 'Name', field: 'name', title: "Name", type: 'String'}
],
rowClasses: [
{className: 'lowlight', field: 'lowlight'}
]}
};
//
// upload functions
//
if (!bidix.upload) bidix.upload = {};
if (!bidix.upload.messages) bidix.upload.messages = {
//from saving
invalidFileError: "The original file '%0' does not appear to be a valid TiddlyWiki",
backupSaved: "Backup saved",
backupFailed: "Failed to upload backup file",
rssSaved: "RSS feed uploaded",
rssFailed: "Failed to upload RSS feed file",
emptySaved: "Empty template uploaded",
emptyFailed: "Failed to upload empty template file",
mainSaved: "Main TiddlyWiki file uploaded",
mainFailed: "Failed to upload main TiddlyWiki file. Your changes have not been saved",
//specific upload
loadOriginalHttpPostError: "Can't get original file",
aboutToSaveOnHttpPost: 'About to upload on %0 ...',
storePhpNotFound: "The store script '%0' was not found."
};
bidix.upload.uploadChanges = function(onlyIfDirty,tiddlers,storeUrl,toFilename,uploadDir,backupDir,username,password)
{
var callback = function(status,uploadParams,original,url,xhr) {
if (!status) {
displayMessage(bidix.upload.messages.loadOriginalHttpPostError);
return;
}
if (bidix.debugMode)
alert(original.substr(0,500)+"\n...");
// Locate the storeArea div's
var posDiv = locateStoreArea(original);
if((posDiv[0] == -1) || (posDiv[1] == -1)) {
alert(config.messages.invalidFileError.format([localPath]));
return;
}
bidix.upload.uploadRss(uploadParams,original,posDiv);
};
if(onlyIfDirty && !store.isDirty())
return;
clearMessage();
// save on localdisk ?
if (document.location.toString().substr(0,4) == "file") {
var path = document.location.toString();
var localPath = getLocalPath(path);
saveChanges();
}
// get original
var uploadParams = new Array(storeUrl,toFilename,uploadDir,backupDir,username,password);
var originalPath = document.location.toString();
// If url is a directory : add index.html
if (originalPath.charAt(originalPath.length-1) == "/")
originalPath = originalPath + "index.html";
var dest = config.macros.upload.destFile(storeUrl,toFilename,uploadDir);
var log = new bidix.UploadLog();
log.startUpload(storeUrl, dest, uploadDir, backupDir);
displayMessage(bidix.upload.messages.aboutToSaveOnHttpPost.format([dest]));
if (bidix.debugMode)
alert("about to execute Http - GET on "+originalPath);
var r = doHttp("GET",originalPath,null,null,username,password,callback,uploadParams,null);
if (typeof r == "string")
displayMessage(r);
return r;
};
bidix.upload.uploadRss = function(uploadParams,original,posDiv)
{
var callback = function(status,params,responseText,url,xhr) {
if(status) {
var destfile = responseText.substring(responseText.indexOf("destfile:")+9,responseText.indexOf("\n", responseText.indexOf("destfile:")));
displayMessage(bidix.upload.messages.rssSaved,bidix.dirname(url)+'/'+destfile);
bidix.upload.uploadMain(params[0],params[1],params[2]);
} else {
displayMessage(bidix.upload.messages.rssFailed);
}
};
// do uploadRss
if(config.options.chkGenerateAnRssFeed) {
var rssPath = uploadParams[1].substr(0,uploadParams[1].lastIndexOf(".")) + ".xml";
var rssUploadParams = new Array(uploadParams[0],rssPath,uploadParams[2],'',uploadParams[4],uploadParams[5]);
var rssString = generateRss();
// no UnicodeToUTF8 conversion needed when location is "file" !!!
if (document.location.toString().substr(0,4) != "file")
rssString = convertUnicodeToUTF8(rssString);
bidix.upload.httpUpload(rssUploadParams,rssString,callback,Array(uploadParams,original,posDiv));
} else {
bidix.upload.uploadMain(uploadParams,original,posDiv);
}
};
bidix.upload.uploadMain = function(uploadParams,original,posDiv)
{
var callback = function(status,params,responseText,url,xhr) {
var log = new bidix.UploadLog();
if(status) {
// if backupDir specified
if ((params[3]) && (responseText.indexOf("backupfile:") > -1)) {
var backupfile = responseText.substring(responseText.indexOf("backupfile:")+11,responseText.indexOf("\n", responseText.indexOf("backupfile:")));
displayMessage(bidix.upload.messages.backupSaved,bidix.dirname(url)+'/'+backupfile);
}
var destfile = responseText.substring(responseText.indexOf("destfile:")+9,responseText.indexOf("\n", responseText.indexOf("destfile:")));
displayMessage(bidix.upload.messages.mainSaved,bidix.dirname(url)+'/'+destfile);
store.setDirty(false);
log.endUpload("ok");
} else {
alert(bidix.upload.messages.mainFailed);
displayMessage(bidix.upload.messages.mainFailed);
log.endUpload("failed");
}
};
// do uploadMain
var revised = bidix.upload.updateOriginal(original,posDiv);
bidix.upload.httpUpload(uploadParams,revised,callback,uploadParams);
};
bidix.upload.httpUpload = function(uploadParams,data,callback,params)
{
var localCallback = function(status,params,responseText,url,xhr) {
url = (url.indexOf("nocache=") < 0 ? url : url.substring(0,url.indexOf("nocache=")-1));
if (xhr.status == 404)
alert(bidix.upload.messages.storePhpNotFound.format([url]));
if ((bidix.debugMode) || (responseText.indexOf("Debug mode") >= 0 )) {
alert(responseText);
if (responseText.indexOf("Debug mode") >= 0 )
responseText = responseText.substring(responseText.indexOf("\n\n")+2);
} else if (responseText.charAt(0) != '0')
alert(responseText);
if (responseText.charAt(0) != '0')
status = null;
callback(status,params,responseText,url,xhr);
};
// do httpUpload
var boundary = "---------------------------"+"AaB03x";
var uploadFormName = "UploadPlugin";
// compose headers data
var sheader = "";
sheader += "--" + boundary + "\r\nContent-disposition: form-data; name=\"";
sheader += uploadFormName +"\"\r\n\r\n";
sheader += "backupDir="+uploadParams[3] +
";user=" + uploadParams[4] +
";password=" + uploadParams[5] +
";uploaddir=" + uploadParams[2];
if (bidix.debugMode)
sheader += ";debug=1";
sheader += ";;\r\n";
sheader += "\r\n" + "--" + boundary + "\r\n";
sheader += "Content-disposition: form-data; name=\"userfile\"; filename=\""+uploadParams[1]+"\"\r\n";
sheader += "Content-Type: text/html;charset=UTF-8" + "\r\n";
sheader += "Content-Length: " + data.length + "\r\n\r\n";
// compose trailer data
var strailer = new String();
strailer = "\r\n--" + boundary + "--\r\n";
data = sheader + data + strailer;
if (bidix.debugMode) alert("about to execute Http - POST on "+uploadParams[0]+"\n with \n"+data.substr(0,500)+ " ... ");
var r = doHttp("POST",uploadParams[0],data,"multipart/form-data; ;charset=UTF-8; boundary="+boundary,uploadParams[4],uploadParams[5],localCallback,params,null);
if (typeof r == "string")
displayMessage(r);
return r;
};
// same as Saving's updateOriginal but without convertUnicodeToUTF8 calls
bidix.upload.updateOriginal = function(original, posDiv)
{
if (!posDiv)
posDiv = locateStoreArea(original);
if((posDiv[0] == -1) || (posDiv[1] == -1)) {
alert(config.messages.invalidFileError.format([localPath]));
return;
}
var revised = original.substr(0,posDiv[0] + startSaveArea.length) + "\n" +
store.allTiddlersAsHtml() + "\n" +
original.substr(posDiv[1]);
var newSiteTitle = getPageTitle().htmlEncode();
revised = revised.replaceChunk("<title"+">","</title"+">"," " + newSiteTitle + " ");
revised = updateMarkupBlock(revised,"PRE-HEAD","MarkupPreHead");
revised = updateMarkupBlock(revised,"POST-HEAD","MarkupPostHead");
revised = updateMarkupBlock(revised,"PRE-BODY","MarkupPreBody");
revised = updateMarkupBlock(revised,"POST-SCRIPT","MarkupPostBody");
return revised;
};
//
// UploadLog
//
// config.options.chkUploadLog :
// false : no logging
// true : logging
// config.options.txtUploadLogMaxLine :
// -1 : no limit
// 0 : no Log lines but UploadLog is still in place
// n : the last n lines are only kept
// NaN : no limit (-1)
bidix.UploadLog = function() {
if (!config.options.chkUploadLog)
return; // this.tiddler = null
this.tiddler = store.getTiddler("UploadLog");
if (!this.tiddler) {
this.tiddler = new Tiddler();
this.tiddler.title = "UploadLog";
this.tiddler.text = "| !date | !user | !location | !storeUrl | !uploadDir | !toFilename | !backupdir | !origin |";
this.tiddler.created = new Date();
this.tiddler.modifier = config.options.txtUserName;
this.tiddler.modified = new Date();
store.addTiddler(this.tiddler);
}
return this;
};
bidix.UploadLog.prototype.addText = function(text) {
if (!this.tiddler)
return;
// retrieve maxLine when we need it
var maxLine = parseInt(config.options.txtUploadLogMaxLine,10);
if (isNaN(maxLine))
maxLine = -1;
// add text
if (maxLine != 0)
this.tiddler.text = this.tiddler.text + text;
// Trunck to maxLine
if (maxLine >= 0) {
var textArray = this.tiddler.text.split('\n');
if (textArray.length > maxLine + 1)
textArray.splice(1,textArray.length-1-maxLine);
this.tiddler.text = textArray.join('\n');
}
// update tiddler fields
this.tiddler.modifier = config.options.txtUserName;
this.tiddler.modified = new Date();
store.addTiddler(this.tiddler);
// refresh and notifiy for immediate update
story.refreshTiddler(this.tiddler.title);
store.notify(this.tiddler.title, true);
};
bidix.UploadLog.prototype.startUpload = function(storeUrl, toFilename, uploadDir, backupDir) {
if (!this.tiddler)
return;
var now = new Date();
var text = "\n| ";
var filename = bidix.basename(document.location.toString());
if (!filename) filename = '/';
text += now.formatString("0DD/0MM/YYYY 0hh:0mm:0ss") +" | ";
text += config.options.txtUserName + " | ";
text += "[["+filename+"|"+location + "]] |";
text += " [[" + bidix.basename(storeUrl) + "|" + storeUrl + "]] | ";
text += uploadDir + " | ";
text += "[[" + bidix.basename(toFilename) + " | " +toFilename + "]] | ";
text += backupDir + " |";
this.addText(text);
};
bidix.UploadLog.prototype.endUpload = function(status) {
if (!this.tiddler)
return;
this.addText(" "+status+" |");
};
//
// Utilities
//
bidix.checkPlugin = function(plugin, major, minor, revision) {
var ext = version.extensions[plugin];
if (!
(ext &&
((ext.major > major) ||
((ext.major == major) && (ext.minor > minor)) ||
((ext.major == major) && (ext.minor == minor) && (ext.revision >= revision))))) {
// write error in PluginManager
if (pluginInfo)
pluginInfo.log.push("Requires " + plugin + " " + major + "." + minor + "." + revision);
eval(plugin); // generate an error : "Error: ReferenceError: xxxx is not defined"
}
};
bidix.dirname = function(filePath) {
if (!filePath)
return;
var lastpos;
if ((lastpos = filePath.lastIndexOf("/")) != -1) {
return filePath.substring(0, lastpos);
} else {
return filePath.substring(0, filePath.lastIndexOf("\\"));
}
};
bidix.basename = function(filePath) {
if (!filePath)
return;
var lastpos;
if ((lastpos = filePath.lastIndexOf("#")) != -1)
filePath = filePath.substring(0, lastpos);
if ((lastpos = filePath.lastIndexOf("/")) != -1) {
return filePath.substring(lastpos + 1);
} else
return filePath.substring(filePath.lastIndexOf("\\")+1);
};
bidix.initOption = function(name,value) {
if (!config.options[name])
config.options[name] = value;
};
//
// Initializations
//
// require PasswordOptionPlugin 1.0.1 or better
bidix.checkPlugin("PasswordOptionPlugin", 1, 0, 1);
// styleSheet
setStylesheet('.txtUploadStoreUrl, .txtUploadBackupDir, .txtUploadDir {width: 22em;}',"uploadPluginStyles");
//optionsDesc
merge(config.optionsDesc,{
txtUploadStoreUrl: "Url of the UploadService script (default: store.php)",
txtUploadFilename: "Filename of the uploaded file (default: in index.html)",
txtUploadDir: "Relative Directory where to store the file (default: . (downloadService directory))",
txtUploadBackupDir: "Relative Directory where to backup the file. If empty no backup. (default: ''(empty))",
txtUploadUserName: "Upload Username",
pasUploadPassword: "Upload Password",
chkUploadLog: "do Logging in UploadLog (default: true)",
txtUploadLogMaxLine: "Maximum of lines in UploadLog (default: 10)"
});
// Options Initializations
bidix.initOption('txtUploadStoreUrl','');
bidix.initOption('txtUploadFilename','');
bidix.initOption('txtUploadDir','');
bidix.initOption('txtUploadBackupDir','');
bidix.initOption('txtUploadUserName','');
bidix.initOption('pasUploadPassword','');
bidix.initOption('chkUploadLog',true);
bidix.initOption('txtUploadLogMaxLine','10');
// Backstage
merge(config.tasks,{
uploadOptions: {text: "upload", tooltip: "Change UploadOptions and Upload", content: '<<uploadOptions>>'}
});
config.backstageTasks.push("uploadOptions");
//}}}