0
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ No, you don't have to use Explicit Euler: DSolve[{x''[t] + 9 x[t] == Cos[3 t], x[0] == 1, x'[0] == 1}, x[t], t] $\endgroup$ – Dr. belisarius Sep 18 '13 at 0:34
  • $\begingroup$ Yes, I know it, but I want compare the exact solution with the approximate solution. But I don't know, do you can help me? please $\endgroup$ – matem Sep 18 '13 at 0:36
1
$\begingroup$

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.}]
$\endgroup$
1
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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