Heun Method - Mathematica Implementation Part 1

Numerical Methods for Solving Differential Equations

Heun's Method

Mathematica Implementation

(continued from last page...)

For reference, I'll repeat the program on this page too:

euler[f_,{x_,x0_,xn_},{y_,y0_},steps_]:=

Block[{ xold=x0,
    yold=y0,
    sollist={{x0,y0}},
    x,y,h
  },

  h=N[(xn-x0)/steps];

  Do[ xnew=xold+h;
    ynew=yold+h*(f/.{x->xold,y->yold});
    sollist=Append[sollist,{xnew,ynew}];
    xold=xnew;
    yold=ynew,
    steps
  ];

  Return[sollist]
]

The red-highlighted parts of the code are the only locations in the program that actually require any modification. If you sift through the rest of the code carefully you should be able to see that it has nothing to do with the algorithm being used, and is merely concerned with the nitty-gritty of constructing our list of solution points, etc., as discussed in great detail in the previous lab.

Obviously the name of the program should no longer be "euler," since we're now planning on using the Heun algorithm. A more appropriate name would be "heun," of course, so we'll make that modification.

More importantly, the method we're using to predict the next y-value is now much more sophisticated, namely:

yn+1 = yn + (h/2)   (f(xn, yn) + f(xn + h, yn +  h f(xn, yn)))

We need to modify the line:

ynew=yold+h*(f/.{x->xold,y->yold});

accordingly, making appropriate changes to put the iteration formula we found for Heun's method into the language of Mathematica. This is something I'm going to leave you to do by yourself. (I warned you that by the end of the course I'd be expecting more independent thinking from you, and this is really not asking too much!) There is actually more than one way to modify the code and get the correct results from the program, some of which even utilize inserting more than one line of new code!

Switch back to Mathematica now, and make the two necessary changes to the Euler program, then save your notebook and come back here to your web browser.

Welcome back. In order to see whether or not your program is working properly, we will now move on and give it a work out.

 

Compass If you're lost, impatient, want an overview of this laboratory assignment, or maybe even all three, you can click on the compass button on the left to go to the table of contents for this laboratory assignment.