Homework 4, due September 22, 2004.

[goodwine]

Unless otherwise indicated, each problem that requires you to write a computer program must be done in C, C++, FORTRAN or another language explicitly approved by the instructor.

1. Consider
   • 
   This is an example of undamped forced vibrations where the forcing frequency and the natural frequency are the same. Here m=1, k=4, F=4 and the forcing frequency is 2.
   1. Determine the analytical solution.
   2. Convert this second order equation into two first order equations.
   3. Write a computer program that uses Euler's method to compute an approximate solution for the time interval t=0 to t=10. Determine a good time step by comparing the approximate solution to the analytical solution.
2. How to use Matlab.
   Matlab has a built-in function called ode45() that is an implementation of the 4th order R-K method. As an example of how to use it, let's use it to numerically solve
   • 
   which is equivalent to the two first order equations
   • 
   To use ode45() you must create a new Matlab function. The easiest way to make a new function is to make a *.m file where the * is the the new function name.
   
   Since we like to solve differential equations with f(x,t) in the equation, let's created a new function called f that needs an x and t as arguments. In this example, the following would go inside the file called "f.m" and the file f.m should be in the same directory in which you started Matlab:

   Code: Select all

   function xdot = f(t,x)
    xdot = zeros(2,1);
    xdot(1) = x(2);
    xdot(2) = sin(t) - x(1) - 0.5*x(2);
   

   Note that the order in the first line, f(t,x) is important. You'll get the wrong answer with f(x,t) there. Now, to compute the solution, simply type the follwing at the Matlab prompt:

   Code: Select all

   >> [t,x] = ode45(@f,[0 10],[1 2])

   This will give back a vector of times and a matrix for x. The first column of x will be x(1) and the second column will be x(2). The [0,10] is the time interval (t=0 to t=10) and the [1 2] are the initial conditions. To plot the answer, just type

   Code: Select all

   >> plot(t,x(:,1))

   Note: if you want to call your function something other than f you have to change f in three places: the first line of the file f.m, the name of the file and the first argument of ode45().
   1. Reproduce the above steps to use ode45() to determine an approximate solution to the given equation (the equation given at the beginning of this problem 2).
   2. Use ode45() to compute a numerical approximation to the solution in problem 1 and compare it with your answer.
3. Consider the famous Lorenz equations
   • 
   where
   • 
   1. Using 4th order Runge-Kutta determine an approximate numerical solution to these equations for the time range t=0 to t=50. Submit a 3D plot of (x,y,z). Be sure to experiment with the time step to ensure that your solution is accurate. The Matlab plot3() function will probably be useful.
      
      Hint: for one first order equation, 4th order R-K uses k1, k2, k3 and k4. To do this for three equations, you will need something like k1, k2, k3, k4, l1, l2, l3, l4, m1, m2, m3 and m4, where the k's are for the first equation, the l's are for the second and the m's are for the third. Furthermore, since the equations for k2, l2 and m2 may require k1, l1 and m1, you will want to compute them in the order of k1, l1, m1, k2, l2, m2, k3, l3, m3, k4, l4, m4.
   2. Use Matlab's ode45() to compute an approximate solution and plot the result.
   3. (5 points extra credit) Explain the significance of these equations.
4. Consider the rather innocent-looking differential equation
   • 
   1. Show that
      • 
      is the solution to the differential equation.
   2. Write a computer program to determine an approximate solution using Euler's method to this first order differential equation.
   3. Submit a plot of the approximate solution and the exact solution for this equation for the time interval from t=-1 to t=1. Use trial and error to determine an appropriate step size by comparing your approximate numerical solution to the exact solution for different step sizes and ensure that the magnitude of the maximum error is less than 0.001. Hint the step size has to be pretty small! Be sure to submit your code as well.
   4. Use the Matlab ode45() function to solve the equation. Submit your code and a plot of the solution and the exact solution. Does Matlab give the correct answer?
   5. (5 points extra credit) Figure out how to tweak the tolerances in ode45(). Is it possible to get an accurate solution?
   6. Write a matlab script to implement the same method you used for part (b) above. (This mistakenly referred to "part (a) above" until 7:00, Monday September 20). An example script using Euler's method to solve the equation
      • 
      would look like the following:

      Code: Select all

      >> x(1) = -6;
      >> dt = 0.01;
      >> t = 0:dt:10;
      >> for i=2:length(t)
      x(i)=x(i-1)+sin(t(i-1))*dt;
      end
      >> plot(t,x)
      

      If you had to use the same size time step as you determined was necessary in part (c), determine approximately how long it will take for Matlab to solve it in this manner.

[grc585]

I have too many passwords to remember.

There are too many extra options on this forum.

I would like something that is easier to print. Why span 4 pages instead of having one homework assignment per webpage? I do not like to click to continue to the next webpage, then click print again.

Now that I am done complaining, I have a real question which I address to Professor Goodwine.

Professor, what was the name of the book (vibrations?) you mentioned in class that is optional, but not required? I would like to get it.

[goodwine]

I have too many passwords to remember.

There are too many extra options on this forum.

I would like something that is easier to print. Why span 4 pages instead of having one homework assignment per webpage? I do not like to click to continue to the next webpage, then click print again.

Then don't print it...

Now that I am done complaining, I have a real question which I address to Professor Goodwine.

Professor, what was the name of the book (vibrations?) you mentioned in class that is optional, but not required? I would like to get it.

There is a link to it on the syllabus. It's optional, meaning it's not required.

[lisaturtle]

I am unsure as to how the code should reflect the Runge-Kutta method in Matlab. I found all the k's, l's, and m's but I don't understand how to find things like x(t+dt/2) so that the program will be able to compute values. I think I'm missing something big here but I don't know what.

[goodwine]

I am unsure as to how the code should reflect the Runge-Kutta method in Matlab. I found all the k's, l's, and m's but I don't understand how to find things like x(t+dt/2) so that the program will be able to compute values. I think I'm missing something big here but I don't know what.

I'm not quite sure which problem you are asking about.

ode45() does it all for you, so you don't have to do any of the k's or anything like that for problems 2 and 3. All you effectively need to do is enter the right hand side of the equation into the function you write. If you are asking about 3 (a), it should be a C or FORTRAN program. I mentioned plot3() to give a hint to how to plot it in Matlab.

For the last problem, you can just use Euler. If you are asking about 4 (f), it should have said the same method as part (b) above instead of part (a).

[mightyduck]

For problem 3b, do you want us to include the number matrices that the matlab code outputs? Because I'm going from 0 to 50 for t, it seems that there will be a lot of pages to print of just numbers.

[mjesick]

For number four, I get x(t) = 1 / ( 1 - e^(-100*t) ) ... a minus sign in the denominator instead of a plus sign

[mjesick]

Because of the absolute value signs in the solution, apparently there are actually four different solutions

[goodwine]

For number four, I get x(t) = 1 / ( 1 - e^(-100*t) ) ... a minus sign in the denominator instead of a plus sign

I'm pretty sure the given solution is correct.

I now see a pontial bit of understandable confusion regarding how I worded the problem. Generally when I (and others) say "show that xxx is a solution" to something, it's o.k. just to substitute it into the equation to show that it satisfies it.

If I say, "determine the solution" or "solve the equation" whether or not I give the answer you need to go through the steps to solve it. So here you can just substitute it to show it satisfies the equation and initial condition.

[goodwine]

For problem 3b, do you want us to include the number matrices that the matlab code outputs? Because I'm going from 0 to 50 for t, it seems that there will be a lot of pages to print of just numbers.

No, just submit a plot. Unless I indicate otherwise, that is sufficient because it shows exactly the same information in a much more efficient manner.

[goodwine]

Some people have reported problems with running over their meager disk quota allowed by OIT. Depending on your quota, problem 4 in this homework may exaserbate this problem since using a very small time step may result in a huge output file.

How to reduce the size of the data file:

1. If you step size is something like 0.00000001, then printing the data to the file every single time step is pointless since you will never see the difference in a plot. The solution is to only plot, say, every 100th time step, or every nth time step, which will reduce the file size by a factor of 100 or n respectively.
2. If you are looking at a plot on the screen or printing a plot, it's also pointless to print the data to a bunch of significant digits. So, you can modify the formatting in the fprintf() function to only print, say, 3 digits after the decimal point. This will also reduce the size of the file by about half since by default it prints, if I remember correctly, 6 digits after the decimal.

Anyway, I modified the following code (Euler integration example from earlier in the course) to do both.

Code: Select all

#include<stdio.h>
#include<math.h>

main() {

  double t,dt,t_start,t_final;
  double x;
  int count=0;
  FILE *fp;

  dt = 0.001;
  t_start = 0.0;
  t_final = 10.0;
  
  x = -1.0;
  
  fp = fopen("data.d","w");

  for(t=t_start;t<=t_final;t+=dt) {
    if(count>=100) {
      fprintf(fp,"%.3f\t%.3f\n",t,x);
      count=0;
    }
    x += sin(t)*dt;
    count++;
  }
  fclose(fp);

  fp = fopen("exact.d","w");
  for(t=0;t<=t_final;t+=0.0001) {
    fprintf(fp,"%.3f\t%.3f\n", t,-cos(t));
  }
  fclose(fp);
}

For what it's worth, this is the sort of problem I'm more than happy to help solve, but expect that later in the course most of you will have developed enough computer experience to handle on your own.

[mightyduck]

I know that although the book had it right, I used the equations for the 4th order R-K from my notebook.... which we found out today were given to us incorrectly a few weeks ago. Any way we could have some points back on our last homework? I'm sure I'm not the only person who made this mistake.

[goodwine]

I know that although the book had it right, I used the equations for the 4th order R-K from my notebook.... which we found out today were given to us incorrectly a few weeks ago. Any way we could have some points back on our last homework? I'm sure I'm not the only person who made this mistake.

What I presented in class was correct, both a few weeks ago as well as today (which is what I also said I was pretty sure was the case in class today). The dt's were in different places, but both are right. Perhaps you homework was incorrectly graded, but this wouldn't be the reason. If you think there was a grading mistake, let me know.

[NDChevy07]

My results for problem 3 quickly blow up to infinity. Is that right, or did I do something wrong?

[goodwine]

My results for problem 3 quickly blow up to infinity. Is that right, or did I do something wrong?

No, it shouldn't blow up for the given parameter values and initial conditions.

[iluvpnutbtr]

I was wonering if we are to plot x1+x2 and the exact or x1 and x2 seperately or just x2 (for problem one)?

[goodwine]

I was wonering if we are to plot x1+x2 and the exact or x1 and x2 seperately or just x2 (for problem one)?

For problem 1, just plot t versus x1. Since x1=x, that's what we care about. Since x2=dx/dt, that plot would be a plot of the velocity versus time. You can plot that if you want, but it isn't necessary.

[iluvpnutbtr]

Help, matlab will not allow me to print or manipulate(insert legend) to any of my graphs. Not just on one computer but on any computer i sign on in the cluseter. i know my plot is corect but i cannot print it out, i just get hunderds of lines of error codes. what is happening?

[goodwine]

Help, matlab will not allow me to print or manipulate(insert legend) to any of my graphs. Not just on one computer but on any computer i sign on in the cluseter. i know my plot is corect but i cannot print it out, i just get hunderds of lines of error codes. what is happening?

I can only guess at the solution to this. Some people have similar problems when the matlab configuration files get corrupted between windows and linux. If you haven't customized any matlab stuff, then one thing that works to solve this is to delete all of your personal matlab config files. To do this run

Code: Select all

rm -rf ~/.matlab/

The config files are hidden in a directory called .matlab/ "rm" is for "remove" and -rf are for "recursive" and "forced" which means go down into any directories and delete it all.

Be careful with the rm command. Until you are used to unix, it can be a bit easier than is advisable to delete all of your files!