本次要求用MATLAB实现二阶微分方程和DNA分析仿真的实验

1. The Pendulum Problem In class, we discussed how the second-order differential equation describing the motion of a simple pendulum can be treated as the two separate differential equations shown

where g is the acceleration due to gravity, 9.81 m/s², L is the length of the pendulum, and θ is the an- gular position relative to the vertical axis. Recall that for a simple pendulum, we neglect complicating factors like the weight of the connecting rod, friction at the hinge, and air resistance. Your task is to write a script to study the motion of such a pendulum.

• Before you start programming, write down the discretized governing equations for this problem using the explicit (forward) Euler method. Use subscripts k and k + 1 to denote the known values at step k and the values to be determined at step k + 1. Remember, when using explicit Euler, the terms on the right-hand side should always reflect information from time step k.
• For a pendulum with length L = 1 m starting from rest position θ₀ = π/3 at t = 0, calculate the angular position, θ, angular velocity, ω, and angular acceleration, α, up to time t = 20 s with

∆t = 0.005 s. Using the velocity and position values, calculate the total energy per unit mass of the pendulum at each timestep. The total energy per mass, E_total, is the sum of potential and kinetic energy as shown below

where h, the vertical displacement from the rest position, is a function of θ and L.

• Plot your three kinematics vectors—position, velocity, and acceleration vs. time—on the same axis. In a separate figure window, plot the total energy of the pendulum vs.
• Repeat this experiment to generate a second set of plots, but this time, treat Equation 2 implicitly as explained in class. That is, the calculation of the angular position at the next step, θ(k + 1), should be based of the newly-calculated velocity at the next step, ω(k + 1), on the right-hand

side as opposed to the current velocity, ω(k). As a reminder, this leads to a “semi-implicit” approach overall since we’re still treating Equation explicitly.

• Using your plots to support your discussion, does forward Euler conserve energy over time? Does using a smaller timestep size, ∆t < 0.005, change this behavior? Include the correspond- ing plots in your report to back up this Does the semi-implicit method tested in Part

(d) conserve energy? Why or why not?

Note: Although you’re encouraged to experiment with different simulation parameters (e.g., length, initial position, timestep size) to generate interesting results for your report, please turn in the version of your code with the values listed in Part (b) and g  = 9.81 m/s² using explicit Euler only.  Don’t forget to include any necessary plots in your report from this week forward!

2. DNA Analysis. In class, we talked about how DNA is good example of how a huge amount of information can be stored in a In this problem, you’ll use iteration and logical expressions to analyze a segment of human DNA stored as a 1-D array. For the sake of programming convenience, the DNA provided to you on CCLE uses numeric values instead of the usual alphabetic abbreviations (A, C, G, and T) as shown in the key below.

Your task is to write a script to calculate the lengths of every protein-coding segment in the DNA as shown in the example below.

C    T    A    A    T    G    A    C    C    T      A    G    A … 2    4     1     1     4     3      1    2     2     4    1     3     1 …

• Begin by downloading the segment of DNA from CCLE (mat) into your current working directory. Next, load the file into MATLAB using the load function as shown below.

load(‘chr1_sect.mat’);

Experiment with this syntax in the command window. Notice that even though there’s no vari- able assignment, MATLAB automatically creates a variable called DNA which your script now has access to.

• Iterate through the codons (the three-base segments shown above) in your DNA vector until you identify the location of the start codon: ATG (or [1 4 3] after conversion). After finding the start codon, continue iterating until reaching the first stop codon: either TAA, TAG, or TGA. Record the length of this protein-coding segment in a new array and begin looking for the next start We’ll ignore any new start codons that may appear before we first encounter one of the stop codons.

Note: Remember, it’s only meaningful to compare intact codons (sets of three) starting from the first element in our vector (we’ll assume this is the only biologically-correct reading frame). For example, dna(1:3) and dna(4:6) are valid groupings; dna(2:4), which splits up two codons, is not valid.

• After reaching the end of the DNA vector, print the following statistics to the command window in the format shown below (don’t panic, these values are for illustration only):

Total Protein-Coding Segments: 1300 Average Length: 99.11

Maximum Length: 1200

Minimum Length: 9

Where the lengths all reflect the number of bases.

• Using your results from Part (c) , what percentage of this DNA is directly used in the protein- coding process (including the start and stop codons)? What is the most frequently and least frequently used stop codon (we’re only interested in counting the stop codons that terminate a protein-coding segment, not any that may be floating around in between these segments)? In the method outlined above, a generic sequence like ATG. . . ATG. . . TAG would yield the number of bases between the first instance of ATG and the stop codon TAG, since we ignore any new start codons that may appear while searching for a stop codon. In detailed words or pseudocode in the report only, how would you modify your algorithm to allow the starting point of a protein- coding segment to be updated such that cases like this yielded the number of bases between the last or most recent instance of ATG and the stop codon?

Note: Please include the chr1_sect.mat in your .zip file when submitting. You’re welcome to use the built in MATLAB function mean, max, and min if you find them helpful. As with all things programming, it’s much easier to start with a small test case where you know the answer first! Try creating your own DNA vector with the numbers in the above examples and confirm that you recover the correct answer before using using the full sequence stored in chr1_sect.mat.

Using the naming convention presented in the syllabus, submit two separate files to the CCLE course web- site: (1) a .pdf of your written report and (2) a .zip file containing all of the MATLAB files written for the assignment. Remember to use good coding practices by keeping your code organized, choosing suitable variable names, and commenting where applicable. Any of your MATLAB .m files should contain a few comment lines at the top to provide the name of the script, a brief description of the function of the script, and your name and UID.