 Consider
 Determine the analytical solution.
 Convert this second order equation into two first order equations.
 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.
 How to use Matlab.
Matlab has a builtin function called ode45() that is an implementation of the 4th order RK method. As an example of how to use it, let's use it to numerically solve
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: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
function xdot = f(t,x) xdot = zeros(2,1); xdot(1) = x(2); xdot(2) = sin(t)  x(1)  0.5*x(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 typeCode: Select all
>> [t,x] = ode45(@f,[0 10],[1 2])
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().Code: Select all
>> plot(t,x(:,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).
 Use ode45() to compute a numerical approximation to the solution in problem 1 and compare it with your answer.
 Consider the famous Lorenz equations
 Using 4th order RungeKutta 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 RK 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.  Use Matlab's ode45() to compute an approximate solution and plot the result.
 (5 points extra credit) Explain the significance of these equations.
 Using 4th order RungeKutta 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.
 Consider the rather innocentlooking differential equation
 Show that
 Write a computer program to determine an approximate solution using Euler's method to this first order differential equation.
 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.
 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 points extra credit) Figure out how to tweak the tolerances in ode45(). Is it possible to get an accurate solution?
 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
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.
Code: Select all
>> x(1) = 6; >> dt = 0.01; >> t = 0:dt:10; >> for i=2:length(t) x(i)=x(i1)+sin(t(i1))*dt; end >> plot(t,x)
 Show that
Homework 4, due September 22, 2004.

 Site Admin
 Posts: 1596
 Joined: Tue Aug 24, 2004 4:54 pm
 Location: 376 Fitzpatrick
 Contact:
Homework 4, due September 22, 2004.
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.
Last edited by goodwine on Thu Jun 15, 2006 4:41 pm, edited 2 times in total.
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.
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.

 Site Admin
 Posts: 1596
 Joined: Tue Aug 24, 2004 4:54 pm
 Location: 376 Fitzpatrick
 Contact:
Then don't print it...grc585 wrote: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.
There is a link to it on the syllabus. It's optional, meaning it's not required.grc585 wrote: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.
Bill Goodwine, 376 Fitzpatrick
Problem 3
I am unsure as to how the code should reflect the RungeKutta 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.

 Site Admin
 Posts: 1596
 Joined: Tue Aug 24, 2004 4:54 pm
 Location: 376 Fitzpatrick
 Contact:
Re: Problem 3
I'm not quite sure which problem you are asking about.lisaturtle wrote:I am unsure as to how the code should reflect the RungeKutta 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.
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).
Bill Goodwine, 376 Fitzpatrick

 Site Admin
 Posts: 1596
 Joined: Tue Aug 24, 2004 4:54 pm
 Location: 376 Fitzpatrick
 Contact:
Re: HW4#4
I'm pretty sure the given solution is correct.mjesick wrote:For number four, I get x(t) = 1 / ( 1  e^(100*t) ) ... a minus sign in the denominator instead of a plus sign
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.
Bill Goodwine, 376 Fitzpatrick

 Site Admin
 Posts: 1596
 Joined: Tue Aug 24, 2004 4:54 pm
 Location: 376 Fitzpatrick
 Contact:
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.chorgan wrote: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.
Bill Goodwine, 376 Fitzpatrick

 Site Admin
 Posts: 1596
 Joined: Tue Aug 24, 2004 4:54 pm
 Location: 376 Fitzpatrick
 Contact:
Disk space issues
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:
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.
How to reduce the size of the data file:
 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.
 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.
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);
}
Bill Goodwine, 376 Fitzpatrick

 Site Admin
 Posts: 1596
 Joined: Tue Aug 24, 2004 4:54 pm
 Location: 376 Fitzpatrick
 Contact:
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.chorgan wrote:I know that although the book had it right, I used the equations for the 4th order RK 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.
Bill Goodwine, 376 Fitzpatrick

 Site Admin
 Posts: 1596
 Joined: Tue Aug 24, 2004 4:54 pm
 Location: 376 Fitzpatrick
 Contact:
Re: Problem 3
No, it shouldn't blow up for the given parameter values and initial conditions.NDChevy07 wrote:My results for problem 3 quickly blow up to infinity. Is that right, or did I do something wrong?
Bill Goodwine, 376 Fitzpatrick
problem one
I was wonering if we are to plot x1+x2 and the exact or x1 and x2 seperately or just x2 (for problem one)?

 Site Admin
 Posts: 1596
 Joined: Tue Aug 24, 2004 4:54 pm
 Location: 376 Fitzpatrick
 Contact:
Re: 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 wrote:I was wonering if we are to plot x1+x2 and the exact or x1 and x2 seperately or just x2 (for problem one)?
Bill Goodwine, 376 Fitzpatrick
matlab
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?

 Site Admin
 Posts: 1596
 Joined: Tue Aug 24, 2004 4:54 pm
 Location: 376 Fitzpatrick
 Contact:
Re: matlab
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 runiluvpnutbtr wrote: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?
Code: Select all
rm rf ~/.matlab/
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!
Bill Goodwine, 376 Fitzpatrick