Mathematical Principles for Scientific Computing and Visualization  
http://www.farinhansford.com/books/scv
by Gerald Farin and Dianne Hansford


This notebook compares Euler's and Huen's methods for solving an initial value ODE problem.  
Problem: Slope field  y' = f(x,y) = y - 0.4x^2 with initial value (-2, 0.76) using 50 sampled points


We compare the solution using 10 steps and 30 steps.
Compare solution with the solution using Huen's method in Huen.nb.

See Chapter 9: Computing Dynamic Processes for details.

Updated May 2008

(* Input a package for vector fields *)

Needs["VectorFieldPlots`"];

(* Set the output folder to be the one this files lives in *)

SetDirectory[ToFileName[Extract["FileName"/.NotebookInformation[EvaluationNotebook[]],{1},FrontEnd`FileName]]];

(* Create a vector field (f_x, f_y) = (1, y - 0.4x^2) over [-2,2]x[-2,2] *)

vecfield=(VectorFieldPlot[{1,y-0.4` (x^2)},{x,-2,2},{y,-2,2},PlotPoints->20,Frame->True,BaseStyle->{"FontFamily"->"Helvetica","FontSize"->11},ColorOutput->CMYKColor]);

(* VectorFieldPlots`  HeadLength->0, *)

(* Get "exact" solution for y' = y - 0.4x^2 for y[-2] = 0.76 over x in [-2,2] *)

sol2=NDSolve[{y'[x]==y[x]-0.4` (x^2),y[-2]==0.76`},y[x],{x,-2,2}];
(* Load function for plotting *)
sol2=sol2[[1,1,2]] ;


(* Euler's method *)

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]]

(* Huen's method *)

heun[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;ynew1=yold+h (f/.{x->xold,y->yold});oldslope=f/.{x->xold,y->yold};newslope=f/.{x->xnew,y->ynew1};ynew=yold+(1/2) h (oldslope+newslope);sollist=Append[sollist,{xnew,ynew}];xold=xnew;yold=ynew,{steps}];Return[sollist]]

(* Solve Euler for y' = y - 0.4x^2 over x in [-2,2], initial value y = -0.76, and 50 steps *)

eulsol=euler[y-0.4` (x^2),{x,-2,2},{y,0.76`},50];

(* Store plot of Euler solution *)

eulplot=ListPlot[eulsol,PlotStyle->{GrayLevel[0.0], PointSize[0.02`]},Frame->True,BaseStyle->{"FontFamily"->"Helvetica","FontSize"->11},ColorOutput->CMYKColor];

(* Solve Huen's method for the same problem *)

heunsol=heun[y-0.4` (x^2),{x,-2,2},{y,0.76`},50];

(* Store plot of Huen's method *)

heunplot=ListPlot[heunsol,PlotStyle->{GrayLevel[0.0], PointSize[0.02`]},Frame->True,BaseStyle->{"FontFamily"->"Helvetica","FontSize"->11},ColorOutput->CMYKColor];

(* Store plot of exact solution *)

splot2=Plot[sol2,{x,-2,2},PlotStyle->{GrayLevel[0.7],Thickness[0.05]},Frame->True,BaseStyle->{"FontFamily"->"Helvetica","FontSize"->11},ColorOutput->CMYKColor];

(* Store initial value point *)

points=ListPlot[{{-2,0.76`}},PlotStyle->{GrayLevel[0.0],PointSize[0.05`]},BaseStyle->{"FontFamily"->"Helvetica","FontSize"->11},ColorOutput->CMYKColor];

(* Plot Vector field, initial value point, Huen solution, Eurler solution, and exact *)
pic=Show[{vecfield,splot2, points,heunplot,eulplot}, Frame->True,BaseStyle->{"FontFamily"->"Helvetica","FontSize"->11}, AxesStyle->Directive[GrayLevel[0.0],Thickness[0.005]],TicksStyle->Directive[GrayLevel[0.0],Thickness[0.005]], ColorOutput->CMYKColor]


(* Create postscript file *)
Export["compare_ODE.eps",pic];