Error in the result

Question

I want solve the following equation:

$ x''(t) + 9 x(t) = Cos(3t) ,\; t \in [0,10] $

$ x(0)= x'(0) =1$

I have to use the Explicit Euler method and a partition of the interval in $N+1$ subintervals (with $N=60$).

I have done this:

f[x_, {u_, v_}] = -9*u + Cos[3 x]
a = 0;
b = 10;
α = 1;
β = 1;
M = 60;
EulerExplicito[f_, α_, a_, b_, M_] := Module[{x, y, h},
    h = N[(b - a)/M];
    x = {a};
    y = {α};
    Do[AppendTo[y, y[[i]] + h*f[x[[i]], y[[i]]]];
       AppendTo[x, x[[i]] + h], {i, 1, M}];
    Return[Transpose[{x, y}]];
    ]

aproxEuler = EulerExplicito[f, {α, β}, a, b, M]
N[aproxEuler]

I don't know what I have to change, but the result is bad.

Answer 1

Compare the exact solution (exacts)

exacts = DSolve[{x''[t] + 9 x[t] == Cos[3 t], x[0] == x'[0] == 1}, 
  x[t], t] 

as belisarius describes in your comments, to explicit Euler (ees) with the code

ees = NDSolve[{x''[t] + 9 x[t] == Cos[3 t], x[0] == x'[0] == 1}, x[t], t,
  Method -> "ExplicitEuler", "StartingStepSize" -> 1/6]

That should be sufficient short of tweaking the subintervals (I think this is the 61 subintervals you are describing). You can then compare a table of your explicit Euler solutions,

Table[{t, Evaluate[x[t] /. ees[[1]]]}, {t, 0, 10, 1/6.}]

and the exact solutions,

Table[{t, Evaluate[x[t] /. exacts[[1]]]}, {t, 0, 10, 1/6.}]

Answer 2

I am not sure what your code is doing, so I prefer to rewrite in a much simple-minded way (in order to fit into my brain). I imagine you need to experiment with simple numerical schemes, so here's how I would tackle the problem. First, I will write down the system of two first order ODEs. Forgive my choice of names for the variables, I am adapting old code...

y'[x] == yp[x]

yp'[x] == -9 y[x] + Cos[3x]

Then I would define the following function for the right hand sides of the above equations (Do not evaluate the above ones!)

f1[x_, y_, yp_] := yp
f2[x_, y_, yp_] := -9y + Cos[3x]

Now I can define a simple Euler procedure. Since I am to compute the solution and its derivative, I will return both (but if you prefer you can only return the pair of coordinates {xn,yn} by adding a compound statement to compute both and then returning the first one only).

nextPoint[{xn_, yn_, ypn_}] := {
    xn + h, 
    yn + h f1[xn, yn, ypn], 
    ypn+ h f2[xn, yn, ypn]}

That's it. All you need now, is to provide the values for the parameters and embed nextPoint into a NestList. Caution: wrap it inside N[] since you do not want to carry tons of symbolic results between steps. This is just the needed initialization

x0 = 0.; n = 600; xn = 10.;
y0 = 1.; yp0 = 1.;

And here we go with the real beef:

p0 = {x0, y0, yp0};
h = (xn - x0)/n;
pts = NestList[ N[nextPoint], p0, n ];

Now you can extract suitable pairs of points. This is the approximate solution:

sol = pts /. {x_, y_, yp_} -> {x, y};
ListPlot[sol, Frame -> True]

and this is its derivative

solp = pts /. {x_, y_, yp_} -> {x, yp};
ListPlot[solp, Frame -> True]

To compare with the exact solution, we first compute it and then evaluate the error in correspondence with the grid points

solrule = DSolve[{u''[x] + 9 u[x] == Cos[3 x], 
    u[0] == u'[0] == 1}, u, x][[1, 1]];
error = pts /. {x_, y_, yp_} -> {x, u[x] - y /. solrule};
ListPlot[error, Frame -> True]

Problem is, this naive implementation of the Euler method does not behave very well. The error becomes awfully large when the number of steps increases (for a given step size) or the step size is too big.

You can adapt the procedure to return the three lists {x0,x1,...,xn}, {y0,y1,...,yn} and {yp0,yp1,...,ypn}, if you prefer to work with one-dimensional vectors of values.