Homework 8, due November 10, 2004.

[goodwine]

Unless otherwise indicated, all problems are from the course text, Elementary Differential Equations and Bounday Value Problem, by Boyce and DiPrima.

1. Section 10.4, numbers 32 and 34.
2. Section 10.5 numbers 18 and 22.
3. Section 10.6 numbers 11 and 12. Note that bounday conditions for 12 were not covered in class, but it's a simple modification of the solution that is clearly covered in the book.
4. Section 10.7 numbers 1 and 6. Note that the exact situation for number 6 was not covered in class, but is a simple enough modification that is covered in the text.
5. Here is the finite difference code for the heat equation covered in class (and corrected for the u_old[] error).

   Code: Select all

   #include<stdio.h>
   #include<math.h>
   #define GRID 100
   #define LENGTH 20
   
   main() {
     double u[GRID+1],u_old[GRID+1];
     double t,dt,t_start=0.0,t_finish=50.0;
     double dx,limit;
     int i,count=1000000;
     FILE *fp;
   
     dx = (double)LENGTH/(double)GRID;
     dt = 0.001;
     limit = t_finish/dt/50.0;       /* only print 50 of the data points */
   
     /* set the initial temperature profile */
     for(i=0;i<GRID+1;i++)
       u[i] = i*dx;
   
     /* set boundary conditions */
     u[0] = 10.0;
     u[GRID] = 0.0;
   
     fp = fopen("data.d","w");
     for(t=t_start;t<t_finish;t+=dt) {    /* loop through time */
       printf("t = %f\n",t);
   
       for(i=0;i<GRID+1;i++) 
         u_old[i] = u[i];                 /* copy the new temp to the old */
   
       for(i=1;i<GRID;i++) {              /* loop through grid */
         u[i] += dt*(u_old[i-1] - 2.0*u_old[i] + u_old[i+1])/pow(dx,2);
       }
   
       if(count>limit) {                  /* only save some of the data */
         for(i=0;i<GRID+1;i++)
   	fprintf(fp,"%f\t%f\n",i*dx,u[i]);
         fprintf(fp,"\n");
         count = 0;
       } else {
         count++;
       }
     }
   
   }
   
   

   Modify the code to compute an approximate numerical solution to problems 18 (a) in Section 10.5 and problem 11 in Section 10.6 (5 points each). Submit your code for each and a plot of the temperature profile at a sufficient number of times to verify your analytical answer and provide an adequate representation of the nature of the solution to the problems.
   Note: I did not yet get a chance to cover this in class, but it must be the case that alpha^2*dt/dx^2 be less than or equal to 1/2 for this finite difference approch to be stable. You can fiddle with GRID and dt if you want, but the values given in the code should work.
6. Write a computer program that determines an approximate numerical solution for the wave equation. Use it to verify your answer for problem 1 in section 10.7. Submit your code as well as plots of the solution that verify your analytical answer as well as effectively communicate the nature of the solution.
   Hint: while not a simple trivial change, this problem is probably best approached by modifying the code for the heat condition equation. Also, be sure to pick dt and dx so that alpha^2*dt/dx is less than or equal to 1.

[student]

I tried to cut and paste the code in #5 to xemacs and compile it in the terminal, but all I got was
/tmp/ccg7T2FO.o(.text+0x204): In function `main':
: undefined reference to `pow'
collect2: ld returned 1 exit status

Did I do something wrong?

[goodwine]

I tried to cut and paste the code in #5 to xemacs and compile it in the terminal, but all I got was
/tmp/ccg7T2FO.o(.text+0x204): In function `main':
: undefined reference to `pow'
collect2: ld returned 1 exit status

Did I do something wrong?

You forgot "-lm" to link the math library (it was saying it doesn't know what the pow() function is).

gcc file.c -lm

[student]

Do we want it to print the values of t from 0-50 in the terminal with the printf command on line 27?

[lisaturtle]

I think I understand the process of finding the displacement for part a, but my answer isn't matching the answer in the back of the book. I checked my work with Mathematica and I get the same thing. I don't understand how there can be sin(n*pi/4) and sin(3n*pi/4) for the k_n's.

[mightyduck]

I'm having troubles with the Mathematica part of the book problems. I typed in

L = 30;
T1 = 30;
T2 = 0;

c[n_] := 60/(n^3 Pi^3) [2 (1 - Cos[n Pi]) - n^2 Pi^2 (1 + Cos[n Pi])]

u[x_, t_, limit_] := 30 - x + Sum[c[n] Exp[-n^2 Pi^2 t/L^2] Sin[n Pi x/L], {n, 1, limit}]

Do[Plot[u[x, t, 50], {x, 0, 30}, PlotRange -> {{0, 30}, {(-5, 25}}], {t, 0, 5}]


This returns nothing but errors and blank graphs... I'm failing to see where I'm going wrong. Help!!

[goodwine]

Do we want it to print the values of t from 0-50 in the terminal with the printf command on line 27?

If GRID is big or dt is small, it can take a while to run, so I put the printf in there just as a means to monitor progress. It doesn't matter whether or not you have that in your program.

[Abigail]

I'm having troubles with the Mathematica part of the book problems. I typed in

L = 30;
T1 = 30;
T2 = 0;

c[n_] := 60/(n^3 Pi^3) [2 (1 - Cos[n Pi]) - n^2 Pi^2 (1 + Cos[n Pi])]

u[x_, t_, limit_] := 30 - x + Sum[c[n] Exp[-n^2 Pi^2 t/L^2] Sin[n Pi x/L], {n, 1, limit}]

Do[Plot[u[x, t, 50], {x, 0, 30}, PlotRange -> {{0, 30}, {(-5, 25}}], {t, 0, 5}]


This returns nothing but errors and blank graphs... I'm failing to see where I'm going wrong. Help!!

In the last line of your posted Mathematica code (the Do[Plot[...]] line) you have an extra '(' in the PlotRange. Take it out, and I think your code should work.

... your friendly TA,
Abigail

[mightyduck]

Sorry. That extra ')' in there was a typo on my part. My actual code does not have it and it still doesn't work. Aaaaah!

[Abigail]

Hmmmm... I typed your code in, and it works for me.

Only difference I can see is your use of square brackets on the RHS defining c[n_]. It seems you've used them as parentheses. That helps with clarity in traditional notation, but in Mathematica it's a bit of a different story. Square brackets are reserved "function" delimiters.

Change those to parens, and see if that works... if not, stop by my office (283 Hurley) sometime this afternoon and we'll take a look.

... your friendly TA,
Abigail

[goodwine]

I'm having troubles with the Mathematica part of the book problems. I typed in

L = 30;
T1 = 30;
T2 = 0;

c[n_] := 60/(n^3 Pi^3) [2 (1 - Cos[n Pi]) - n^2 Pi^2 (1 + Cos[n Pi])]

u[x_, t_, limit_] := 30 - x + Sum[c[n] Exp[-n^2 Pi^2 t/L^2] Sin[n Pi x/L], {n, 1, limit}]

Do[Plot[u[x, t, 50], {x, 0, 30}, PlotRange -> {{0, 30}, {(-5, 25}}], {t, 0, 5}]


This returns nothing but errors and blank graphs... I'm failing to see where I'm going wrong. Help!!

It took me a while, but I found your error. It's the square brackets in your definition of c[n]. It should be

Code: Select all

c[n_] := 60/(n^3 Pi^3) (2 (1 - Cos[n Pi]) - n^2 Pi^2 (1 + Cos[n Pi]))

In Mathematica, square brackets [] are used for functions, e.g., Sin[t], and regular parentheses are used to group mathematical operations, e.g., 3*(5+6).

Addendum: I didn't notice the above post. It has the same answer -- drop the square brackets except for functions.

[lisaturtle]

Disregard my stupid post...turns out that I was doing the wrong problem, but at least I was getting it right!

[goodwine]

I think I understand the process of finding the displacement for part a, but my answer isn't matching the answer in the back of the book. I checked my work with Mathematica and I get the same thing. I don't understand how there can be sin(n*pi/4) and sin(3n*pi/4) for the k_n's.

Sorry it has taken me so long to reply. I'm definitely getting an answer that I can't correlate to the back of the book. However, I plotted du/dt(x,0) and it was clearly g(x), so the back is o.k.

One thing you can do now, though, is to plot u'(x,0} for your answer and you should get g(x).

[Grendel]

For 10.6.11 , i see that cn changes with even and odd n. Is there a way to tell mathematica to use 0 for odd and -4/(n^2-1)*Pi ??

I tried:

Code: Select all

c = If [ OddQ[n] , 0 , (-4*Pi) / (n^2-1) ]

I tried using this with a Plot3D of the answer given in the back of the book and i got error message "1/0 encountered". is there a way to make this work?

[Grendel]

Nevermind, i got it to work...

[goodwine]

Nevermind, i got it to work...

You're making my job too easy...

[mightyduck]

For problem 6, are we just making a 2nd order equation into two first order equations? If so, I'm still not doing it correctly. Perhaps a little direction, because all I get are very extremely large numbers or that little "NaN" thing.

[goodwine]

For problem 6, are we just making a 2nd order equation into two first order equations? If so, I'm still not doing it correctly. Perhaps a little direction, because all I get are very extremely large numbers or that little "NaN" thing.

Yes, you have the right idea. The only difference between the wave equation and heat conduction equation is that the time derivative is second order in the former. The consequence of this fact in the finite difference numerical approximation will be that you will have a second order differential equation for u.

Let's use the syntax from the provided heat equation code. Let u be the string displacement at location i*dx (you can change it to y if you want). Since we have a second order equation for u, we need to convert it (for each i!) to two first order equations. Let's do u and v (u for displacement and v for velocity).

Thus, at location i*dx, you have something close to

u'' = (u[i+1 - 2*u + u[i-1])/pow(dx,2)

(there's probably an alpha going on in there as well). Since v = u', v' = u''[i], which we know. Thus using Euler, the relevant lines in the code will look something like

v[i] += (u_old[i+1] - 2*u_old[i] + u_old[i-1])*dt;
u[i] += v_old[i]*dt;

(and at each time step you should copy u[i] and v[i] into u_old[i] and v_old[i]).

Summary and quick answer: since it's second order, you need to make a "velocity" array as well as an "old_veolocity" array and for each i, update the string position in two steps (one for position and one for velocity) just like you did for a mass-spring-damper system.

[7]

For 10.5:18, I am running the C code and it seems to be stopping at 29.990 seconds and not 30.0, so I am not getting the right answer in the back of the book. Am I doing something wrong?

[asayers]

v += (u_old[i+1] - 2*u_old + u_old[i-1])*dt;
u += v_old*dt;

Should the v term be divided by (dx)^2?

[goodwine]

v += (u_old[i+1] - 2*u_old + u_old[i-1])*dt;
u += v_old*dt;

Should the v term be divided by (dx)^2?

Yes.

[goodwine]

For 10.5:18, I am running the C code and it seems to be stopping at 29.990 seconds and not 30.0, so I am not getting the right answer in the back of the book. Am I doing something wrong?

In a for() loop, the loop continues as long as what is after the first semicolon is true. It's probably something like t<t_final. In your case, if dt is greater than 0.01, then it won't do the next iteration since that would make t>t_final.

I just eyeballed my plots, but they were very close to the answer in the back for all three values for alpha, so it should be true for the first one as well. I cannot remember what I was using for dt or for GRID when I did it. Perhaps you need to increase the resolution, i.e., increase GRID and decrease dt (make sure the necessary relationship between dt and dx is true, though!).

Your plots should look like you would expect: zero at the two ends and curved to a maximum value in the middle and the overall temperature distribution should decrease with time.

[HugsyMcHugs]

Are the solutions to either HW 7 or HW 8 going to be posted tonight to use for the test?

[Abigail]

Are the solutions to either HW 7 or HW 8 going to be posted tonight to use for the test?

I'll be posting solutions to Homework 7 later this evening... don't touch that dial!

Dr. Goodwine has posted partial solutions to Homework 8. I'll post another part of the solutions later tonight; between us, we may or may not cover the whole assignment.

Complete solutions to Homework 8 unfortunately won't be ready to post until after the exam. Sorry about that...