Mathematica 1, Due 11:59pm, 8/27 on paper

      0.1, 0.9 using Mathematica.

      0.11, 0.15, 0.19 using pen and paper.

    C Programming 1, Due 11:59pm, 9/3

      (0.25) Write a C (or C++) program to sum the series

        1+1/2+1/3+1/4+1/5+1/6+......
      Mathematically the series is divergent. Is your result divergent? If not, why?

    Mathematica 2, Due 11:59pm, 9/3 on paper

      1.3, 1.5, 1.11, 1.29, 1.32, 1.60, 1.74, 1.79
      Using Mathematica.

    C Programming 2, Due 11:59pm, 9/10

      Write a C program to solve the equation, x^2+x-1=0, using both the fixed-point iteration and the Newton's method. Compute the result up to 10 significant digits. The program should be written in such a way that it will take the initial guess as an input from the terminal, and return the root together will the number of iterations taken. If the iteration takes more than 100 times then stop and output "Not convergent after 100 iterations." Compare the result produced by two methods, and show that the order of convergence is 1 and 2 respetively.

    Mathematica 3, Due 11:59pm, 9/10 on paper

      2.1, 2.2, 2.3, 2.4.
      Solve equations in problems 2.5, 2.6.
      2.27, 2.28.

    C Programming 3, Due 11:59pm, 9/17

      The program pdec.c solve the general matrix equation A*X=B using the LU decomposition. Modify the program such that in addition to solution of the equation, it also calculats the det(A), tr(A), and the inverse matrix of A.

    Mathematica 4, Due 9:30am, 9/24 on paper

      2.61.
      First plot the two curves to locate the approximate solution near (x=1.6,y=-2.3), Then you can use either Newton's or the fixed-point iteration method to find the root up to 5 significant digits. (If you are using the fixed-point method, solve the first equation for y and the second equation for x.)

      3.2, 3.34, 3.35, 3.73.
      For problem 3.73, you can use the least-square-fitting function provided by the Mathematica.

      4.14, 4.36, 4.44.
      Voluntary: Go through all exercises in Chapter 4, Section 4.15 (page 372 -376).

    C Programming 4, Due 11:59pm, 10/1

      The radioactive intensity of a given amount of radioactive material follows the exponential decay, I(t)=A*exp(-t/T), where t is the time, A and T are constants. In a lab experiment a set of DATA is observed. Assuming A=7.5 is known, write a C program to find the T using the least-square fit method. The program must be able to read data from the file "data". The first line in "data' is the total number of points which is <= 200.

      The half-life T_0 is defined as I(T_0)=I(0)/2. How many days is the half-life of Radon?

      Extra Credit: Estimate the uncertainty in your answer.

    Mathematica 5, Due 9:30am, 10/1

      1. Using the Euler method to solve equation dy/dx=x/y , y(x=0)=1.0, in the interval (0,1), using the step size h=0.1. Compare your result with exact solution y=sqrt(1+x^2) by plotting both.

      2. Repeat the problem using the 4th order Runge-Kutta method.

      3. Varies the step size h to show that the global error for two methods scales as O(h) and O(h^4). (The global error can be estimated by comparing the numerical solution at y(x=1.0) with the exact result.)

    Project 1: Pendulum Dynamics, Due 10/14
    You can do this project either in C or Mathematica

      Use the 4th order Runge-Kutta method explicitly to solve these problems:

      1. Simple Pendulum: d^2 x(t)/d^2 t = -sin(x)
        Obtain a phase space portrait similar to Figure 12.3 (all figures refer to the handout from De Jong's Book).
      2. Damped Pendulum: d^2 x(t)/d^2 t = -sin(x)- c (d x/ d t)
        Obtain a phase space portrait similar to Figure 12.4.
      3. Damped Driven Pendulum: d^2 x(t)/d^2 t = -sin(x)(1+F cos(w t))- c (d x/ d t)
        Obtain a phase space portrait similar to Figure 12.5. Take w=2, c=0.1, h=pi/20.
      4. The Road To Chaos: Obtain phase space portraits for limit cycle, periodic doubling, and doubling again (Figures 12.6,7,8, the d value given in these figures is F in the text.). (You may need as many as 10,000 steps.)
      5. Chaotic Orbit: Obtain a phase space portrait similar to Figure 12.9.
      6. Pointcare Section: Construct the Pointcare section for the chaotic pendulum (Figure 12.11).
      7. Duffing Oscilattor: Modified your program to solve the Duffing Two-well oscillator in the chaotic regime. Constructe the Pointcare section and demonstrate the self-similarity the of thestrange attractor (Figure 12.12, 12.13, 12.14).

      Note 1: Make sure that your program runs for a loop with only few steps before you let it runs for many step.

      Note 2: Part 5,6,7 may require very long run, C programming is recommended.

    Mathematica/C 6, Due 9:30am, 10/22

      Using the Monte Carlo method to solve the following problem.

      Three unit circles (with radius 1 ) are arranged in such a way that their centers formed a equilateral triangle with the distance betweem the center being d. Calculate their common overlap area as d varies from 0 to 2. The result should be accurate for two significant digits.
      Note 1: show your result as a plot of the area vs d.
      Note 2: Make a brief description on how you solve the problem.

    Mathmatica/C Program 7: Temperature Distribution in a Plate, Due 10/29

    A rectangular plate is defined by four boundaries: y=0, y=1.0, x=0, x=1. At time t=0, the temperature distribution u(x,y) at the boundary is fixed as: u(x=0,y)=300*y, u(x=1,y)=300*y, u(x,y=0)=0, u(x,y=1)=800*(x-0.5)^2+100.

    1. Steady State (Laplacian Equation)
      Using iterative method to solve the steady state temperature distribution. The minimum grid size should be 20X20 (namely x=0,0.1,0.2,...,1.0; y=0,0.1,0.2,...,1.0). The iteration should be carried out such that the results is good up to 2 significant digits. The answer should be shown in a 3-D picture landscape plot and a 2-D isothermal contours plot.

    2. Time Dependence (Parabolic PDE)
      Using the simple difference method to solve the time dependence of the temperature distribution. (Hint: use the formula given in page 636.) Choose r=0.1, k/(c*rho)=0.01. Beside the boundary condition given in part I, the initial condition is u(x,y,t=0)=0 at the interior of the plate.
      a) Obtain u(x,y,t=1) in 3-D plot.
      b) Obtain u(x,y,t=100) in a 3-D plot, compare it with part I.
      c) Plot the temperature at the center as a function of time from t=0 to t=100.

    Submit a hard copy of program with your figures. Make sure you label all figures properly and clearly.

    C Program 8: Drum beats -- 2-D Wave Equation, Due 11/10.

    The vibration of a drum surface is governed by the 2-D wave equation. For simplicity let's consider a square drum of dimension 1X1 meter, and take the sound velocity to be c=1 m/s. (Pick any dx you want, but you probably need 40X40 grid to get reasonable result.) The perimeter of the drum is fixed to a steady frame. At time t=0, the drum is hit at the center, with the initial disturbance 

    u(x,y,t=0)=0.1 exp(-200*(x^2+y^2))
    and zero velocity. Investigate the dynamics of drum by solving the 2-D wave equations.

    1. Obtain a 3D plot of the drum surface at following time: t=0.1, 0.25, 0.5, 0.75, 1.0 s.
    2. Discuss qualitatively whether or not your results agree with intuition.
    3. Solve the same problem with circular boundary condition. Namely, the drum surface is fixed at the boundary defined by the equation x^2+y^2=0.5^2. Discuss the difference between the two boundary conditions.
    Project 2: Computerized Tomography, Due 11/24.

    This is your last assigntment, but also the most difficult one. After all, we are dealing with a real life problem. It counts for 15% of your final grade. (Extra credit will be given for exceptional performance). You are encouraged to discuss with each others, and with Brian and me.

    Computerized Tomography (CT) is one of the most potent diagnosis tools in modern medical technology. In CT a seriers of projection funcion ct_projection(xi,phi) is measured, where the projection direction is perpendicular to the axis xi and phi is the angle between the xi axis and the object x axis. The basic objective of CT is to reconstruct the image from these projection functions. (Note: the pojection functions ct_projection(xi,phi), unknown_projection(xi,phi), circle_projection(xi,phi), dounut_projection(xi,phi) and the FFT routines fft61(double data[ ],int n,int sign) are all contained the program projection.c )

    I. Write a program to do the following:

    1. For a given fixed phi (for example, phi=0), sample the projection function ct_projection(xi,phi) uniformly N (must be a power of 2) time to form the projection array proj_fun(xi_i) for $-1<=xi_i<1$. Plot proj_fun(xi_i)
    2. Use the fast Fourier Transform procedure fft61(proj_fun,N,0) to perform the discret Fourier transform. Plot the transformed projection function proj_fun(k_i) in k space.
    3. Use fft(proj_fun,N,1) to perform the inverse Fourier transform on proj_fun(k_i) to recover proj_fun(xi_i). Compare it with what you start with.

    II. Write a program to perform CT using FFT. Follow the modulus steps listed in the handout:

    1. Construct proj_fun(k_i) from ct_projection(xi,phi) for a given phi.
    2. Perform fft on proj_fun(k_i) to obtain big_proj(k_i).
    3. Perform inverse FFT on big_proj(k_i)*|k| to obtain the modified projection function mproj_fun(xi_i) .
    4. calculatin the contribution to the image object1(x,y) for every pixel (x_i,y_j).
    5. Repeat steps 1-4 for several different projection directions. Compare the accumulated image with the image function ct_objects.dat.

    III. Repeat part II by sampling more projection directions until your image agrees with the original (use your judgement). How many projection directions do you need?

    IV. Given three different projection functions unknown_projection(xi,phi), circle_projection(xi,phi), nounut_projection(xi,phi), use your program to reconstruct the three images.