Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I want to solve this equation with NDSolve[] using the Euler method:

x'[t] == 0.5*x[t]-0.04*(x[t])^2

with initial condition and step size

x[0] == 1, h == 1 

and final time = 10.

I also don't know how to make a table including x'[t] and x[t] and time.

share|improve this question

3 Answers 3

up vote 13 down vote accepted

NDSolve has a slew of options that allow you to control the method. You can find the standard reference here.

There, we learn how to access Euler's method using NDSolve:

Clear[x];
x = x /. First[
  NDSolve[{x'[t] == 0.5*x[t] - 0.04*(x[t])^2, x[0] == 1}, x, {t, 0, 10},
  StartingStepSize -> 1, Method -> {"FixedStep", Method -> "ExplicitEuler"}]
];
grid = Table[{t, x[t]}, {t, 0, 10, 1}]
ListLinePlot[grid]

plot of differential equation solution

It is also quite easy to program this from scratch, if you prefer.

share|improve this answer
6  
No need for "FixedStep"; since Mathematica already "knows" that the Euler method is a fixed step method, one can do this instead: NDSolve[{x'[t] == 0.5*x[t] - 0.04*(x[t])^2, x[0] == 1}, x, {t, 0, 10}, StartingStepSize -> 1, Method -> "ExplicitEuler"] –  J. M. Oct 11 '12 at 23:05

As Mark wrote "It is also quite easy to program this from scratch, if you prefer."

For a solution without NDSolve try..

MyEuler[start_, end_, initialvalue_, nrOfsteps_] :=Module[
       {a = start, b = end, j, m = nrOfsteps},
       h = (b - a)/m; (* fixed step-size *)
       T = Table[ a + (j - 1) h, {j, 1, m + 1}]; 
       Y = Table[ initialvalue, {j, 1, m + 1}];  
       For[ j = 1, j <= m, j++,
            Y[[j + 1]] = Y[[j]] + h  f[T[[j]], Y[[j]]];
           ]; 
       Transpose@{T, Y}] 

Testing

f[t_, x_] = 0.5*x - 0.04*x^2;(* rhs of your ODE *) 
pts = MyEuler[0.0, 10.0, 1.0, 10]; 
ListLinePlot[pts, Mesh -> All, MeshStyle -> Red, Frame -> True]

enter image description here

share|improve this answer
2  
something like myEuler[t0_, t1_, initial_, steps_] := With[{h = (t1 - t0)/steps}, FoldList[#1 + h f[#2, #1] &, initial, Range[t0 + h, t1, h]]] could do –  belisarius Oct 11 '12 at 20:53
    
@belisarius Sure! But the OP wanted some Table. I guess by table he meant the array to store the discrete states as commonly done in most textbooks. I felt that mundane For will also better suit a classic textbook Euler .. –  PlatoManiac Oct 11 '12 at 21:09
    
You could also go one step further and show how to exten NDSolve with your method. –  user21 Oct 11 '12 at 23:24
    
@ruebenko I think J.M already has done exactly that while we were sleeping in Europe..right? –  PlatoManiac Oct 12 '12 at 7:04
    
yes, I think he must have posted this a few minutes after the comment... –  user21 Oct 12 '12 at 10:19

Had the Euler method not been built-in, one could still use NDSolve[]'s method plug-in framework, which enables NDSolve[] to "know" how to use Euler's method.

Here's how to "teach" NDSolve[] the Euler method:

Euler[]["Step"[rhs_, t_, h_, y_, yp_]] := {h, h yp};
Euler[___]["DifferenceOrder"] := 1;
Euler[___]["StepMode"] := Fixed;

Plugging in the "new" method into NDSolve[] is a snap:

xa = x /. First @ NDSolve[{x'[t] == 0.5*x[t] - 0.04*(x[t])^2, x[0] == 1}, 
                          x, {t, 0, 10}, Method -> Euler, StartingStepSize -> 1];

Getting the corresponding table is easily done, thanks to the special methods for accessing the internals of an InterpolatingFunction[]:

pts = Transpose[Append[xa["Coordinates"], xa["ValuesOnGrid"]]]
   {{0., 1.}, {1., 1.46}, {2., 2.10474}, {3., 2.97991}, {4., 4.11467}, {5., 5.49478},
    {6., 7.03447}, {7., 8.57235}, {8., 9.91912}, {9., 10.9431}, {10., 11.6246}}

Showing the InterpolatingFunction[] and the points together in one plot is also easily done:

Plot[xa[t], {t, 0, 10}, Epilog -> {AbsolutePointSize[4], Red, Point[pts]}, Frame -> True]

solution of differential equation via Euler

If you wanted the derivatives as well in your table, it's easy to modify pts:

phs = Append[#, xa'[#[[1]]]] & /@ pts;

You can use these for the "phase plot" of your differential equation:

ParametricPlot[{xa[t], xa'[t]}, {t, 0, 10}, AspectRatio -> 1/GoldenRatio, 
               Epilog -> {AbsolutePointSize[4], Red, Point[Rest /@ phs]}, 
               Frame -> True]

phase plot

The phase plot looks ugly here, since the step size was not very small and the interval of integration was not sampled well. Due to this, the InterpolatingFunction[] is unable to make a smooth-looking derivative. A smaller value of StartingStepSize (say, $1/20$) would have resulted in something that looks smoother:

smoother phase plot

but of course at the expense of more evaluations of the right-hand side.

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.