Homework 8, due Wednesday November 1.

[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. Write a computer program to determine an approximate numerical solution to
   • 
   using Euler's method, 2nd order R-K (called the "improved Euler method" in the course text) and 4th order R-K. Then cut the time step in half. Does the local truncation error after the first time step behave as expected? Does the global accumulated error behave as expected?
2. Write a computer program to determine an approximate numerical solution for 0<t<1.5 to
   • 
   Determine an approrpriate time step by either determining tne analytical solution or by continuing to decrease the time step until it appears that the solution no longer changes.
3. 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().
   • Use ode45() to find approximate numerical solutions for the differential equations in the first two problems and compare the result with your answer. You don't have to plot them on the same plot, but at least visually compare them to make sure you get the same answer.
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 either the 2nd order or 4th order R-K 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 Euler's method to determine an approximate numerical solution to the above "innocent" equation. 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)
      

      Using this script and Euler's method. Approximately determine how long it would take for the script to run to return an approximate solution that is as accurate as your code from above.

[cverhuls]

Whenever I try to graph #4 using 4th order R-K, the value of x never changes (even for dt=.001). I'm guessing it's due to the fact that the slope is zero except for the very sharp increase around t=0.

The ode45 function in MATLAB gives the same results.

Knowing that x(0)=.5 (from the exact solution), I can plot the graph in two parts. However, this uses more than the given diff eq and initial condition.

How should I fix this?

[goodwine]

Whenever I try to graph #4 using 4th order R-K, the value of x never changes (even for dt=.001). I'm guessing it's due to the fact that the slope is zero except for the very sharp increase around t=0.

The ode45 function in MATLAB gives the same results.

Knowing that x(0)=.5 (from the exact solution), I can plot the graph in two parts. However, this uses more than the given diff eq and initial condition.

How should I fix this?

Keep trying smaller and smaller dt's.

Matlab gives the wrong result. For me it wasn't constant, though. It just went from -1 to +1 at too late of a time.

[phicks]

For Problem 1,
I am getting the local truncation error for all methods to be roughly four, and overall the Euler Method is the closest, followed by rk4, which is followed closely by rk2. I have tried a number of different time steps and I used both real and double precision but got the same results. I am fairly certain that my code is correct. I wonder if this equation could be one of those weird ones where you don't get what you expect?

Thanks, Phillip

[goodwine]

For Problem 1,
I am getting the local truncation error for all methods to be roughly four, and overall the Euler Method is the closest, followed by rk4, which is followed closely by rk2. I have tried a number of different time steps and I used both real and double precision but got the same results. I am fairly certain that my code is correct. I wonder if this equation could be one of those weird ones where you don't get what you expect?

Thanks, Phillip

I checked using Euler and 4th order R-K and the local truncation error changes approximately as expected for me.

Are you saying that the actual error after the first step is 4? Since the initial condition is 1, your error is 4 times larger than the magnitude of the solution. You either have an error in your code or your time step is way too big. FWIW, I used .5 and .25 for time steps in my comparison.

[pcourtne]

I am getting the same thing as phicks. My Euler method approximation is the closest to the exact answer and I've used a bunch of different time steps.

[goodwine]

I am getting the same thing as phicks. My Euler method approximation is the closest to the exact answer and I've used a bunch of different time steps.

The local truncation error (the error after one time step) when I use dt=0.5 is 100 times smaller for the 4th order R-K than for Euler when I do it.

[goodwine]

I am getting the same thing as phicks. My Euler method approximation is the closest to the exact answer and I've used a bunch of different time steps.

A couple other considerations:

• Is your exact answer correct?
• Is dt less than 1? If it's larger than 1 then higher order methods could have greater error.
• Since the solution is unstable, i.e., blows up, for how long are you running it? I think up to t=2 or 3 should suffice.

I'm doing Euler and 4th order R-K right now and R-K is much better except for larger dt.

[dlipp]

Prof. Goodwine,

I keep getting the following error when trying to load my data from the Fortran file into MATLAB:

??? Error using ==> load
Number of columns on line 13 of ASCII file /afs/nd.edu/user19/dlipp/Private/Diff Eq/rk4.d must be the same as previous lines.

It's not always line 13 that is the issue. Sometimes it's 1, sometimes 100. I printed the output in the terminal and found it to be when time is equal to zero. I defined the time parameters and then did a format and that's when it jumped to line 13. Any pointers?

[goodwine]

Prof. Goodwine,

I keep getting the following error when trying to load my data from the Fortran file into MATLAB:

??? Error using ==> load
Number of columns on line 13 of ASCII file /afs/nd.edu/user19/dlipp/Private/Diff Eq/rk4.d must be the same as previous lines.

It's not always line 13 that is the issue. Sometimes it's 1, sometimes 100. I printed the output in the terminal and found it to be when time is equal to zero. I defined the time parameters and then did a format and that's when it jumped to line 13. Any pointers?

It's saying that each row in the file does not have the same number of numbers in it. The probable cause for this is to have more than one write() statement in the program, but the different write statements do now write the same number of things.

[cchinske]

Professor Goodwine, did you receive my email regarding the Linux machine error?

Christopher Chinske

[Kalbante]

I got the same ASCII error when only when I had huge data files.
It seems like after a certain size the ASCII file or Matlab makes an error

[goodwine]

I got the same ASCII error when only when I had huge data files.
It seems like after a certain size the ASCII file or Matlab makes an error

If some problems are caused by huge data files, the way to solve that is to put some if-then statements in the program to only print, say, every 1000th data point to the file. If the data file is so huge that it takes a while to load into matlab or causes other errors, then even if you only plot every 1000th (or whatever) point you won't be able to tell.