3
$\begingroup$

According to what I have read, the first formula is the classical Euler method and the second is the improved Euler method for second-order equations.

  • Method A: accuracy of order h
S[a_, b_, h_, N_] := (u[0] = a; u[1] = a + h*b; 
  Do[u[n + 1] = 
    2 u[n] - u[n - 1] + h*h*f[n*h, u[n], (u[n] - u[n - 1])/h], {n, 1, 
    N}])
  • Method B: accuracy of order h^2
Q[a_, b_, h_, N_] := (u[0] = a; v[0] = b; 
  Do[{u[n + 1] = 
     u[n] + h*
       F[u[n] + (h/2)*F[u[n], v[n]], 
        v[n] + (h/2)*
          G[u[n], v[
            n]]],                                                     \
          
    v[n + 1] = 
     v[n] + h*
       G[u[n] + (h/2)*F[u[n], v[n]], v[n] + (h/2)*G[u[n], v[n]]]}, {n,
     0, N}])

All the examples I have seen both in the book I have and in YouTube videos deal with two function systems and initial values. My question is this if we are given more than two functions how should we work? I can't find an example guide to figure out how. Any example would be appreciated. Thank you in advance

For example

It is given the following problem: $$X'=Z+(Y-\alpha)X$$ $$Y'=1-\beta Y-X^2$$ $$Z'=-X-\gamma Z$$ with initial conditions $(X(0),Y(0),Z(0)=(1,2,3)$.

Where

X: interest rate

Υ:investment demand

Z: price index

$\alpha$: savings, $\beta$: cost per investment, $\gamma$: the absolute value of the elasticity of demand

And we want the results with the Improved Eulers method. It is obvious that I have to use Method B: accuracy of order h^2. But I do not know how to define the new function on Mathematica

Improve this question
$\endgroup$
8
  • 2
    $\begingroup$ To me, this sounds more of a Mathematics Question, appropriate for Math SE. $\endgroup$
    – Domen
    Commented Jan 26, 2023 at 23:46
  • $\begingroup$ @Domen Fyi I have updated my post so you will understand what I mean $\endgroup$
    – Athanasios Paraskevopoulos
    Commented Jan 26, 2023 at 23:56
  • $\begingroup$ @cvgmt I created a new question about Improved Eulers Method as we discussed earlier today. $\endgroup$
    – Athanasios Paraskevopoulos
    Commented Jan 27, 2023 at 0:27
  • $\begingroup$ 1. 1st implementation you've found is bounded on second order equations i.e. equation in form $y′′=f(x,y,y′)$ so it cannot be extended to first order system i.e. $y′=f(x,y)$, at least cannot be extended in a straightforward way. 2. Euler's method can be defined based on first order system in a much simpler manner, and n-th order system can always be transformed to 1st order system, see e.g. mathematica.stackexchange.com/a/158519/1871, notice 2nd implementation you've found is essentially based on this idea. BTW… $\endgroup$
    – xzczd ♦
    Commented Jan 27, 2023 at 2:59
  • $\begingroup$ …it's usually called midpoint method: en.wikipedia.org/wiki/Midpoint_method 3. 2nd implementation you've found isn't convenient for extending, to implement it in a more convenient way, observe the output of Sin[{Pi/3, Pi/6}] and think about how you can make use of this feature of Mathematica. 4. Please first make some effort to understand Euler's method itself, the wiki page isn't a bad source: en.wikipedia.org/wiki/Euler_method $\endgroup$
    – xzczd ♦
    Commented Jan 27, 2023 at 3:07

1 Answer 1

Reset to default
3
$\begingroup$

This works but it is not as accurate as the NDSolve result

α=0.9;
β=0.2;
γ=1.2;

f[x_,y_,z_]:=z+(y-α)*x
g[x_,y_,z_]:=1-β*y-x^2
p[x_,y_,z_]:=-x-γ*z

Q[a_,b_,c_,h_,N_]:=(u[0]=a;v[0]=b;w[0]=c;
  Do[{u[n+1]=u[n]+h*f[u[n]+h/2*f[u[n],v[n],w[n]],
                      v[n]+h/2*g[u[n],v[n],w[n]],
                      w[n]+h/2*p[u[n],v[n],w[n]]],
      v[n+1]=v[n]+h*g[u[n]+h/2*f[u[n],v[n],w[n]],
                      v[n]+h/2*g[u[n],v[n],w[n]],
                      w[n]+h/2*p[u[n],v[n],w[n]]],
      w[n+1]=w[n]+h*p[u[n]+h/2*f[u[n],v[n],w[n]],
                      v[n]+h/2*g[u[n],v[n],w[n]],
                      w[n]+h/2*p[u[n],v[n],w[n]]]},
      {n,0,N}]);

Q[1,3,2,0.1,1000]

X=Interpolation[Table[{n,u[n]},{n,0,1000}]];

Y=Interpolation[Table[{n,v[n]},{n,0,1000}]];

Z=Interpolation[Table[{n,w[n]},{n,0,1000}]];

ParametricPlot3D[{X[t],Y[t],Z[t]},{t,0,1000}]

enter image description here

Improve this answer
$\endgroup$
1
  • 1
    $\begingroup$ Not too bad in my view. You can improve it by making the grid denser, try e.g. h = 0.01; tend = 100; nmax = tend/h; Q[1, 3, 2, h, nmax] // AbsoluteTiming $\endgroup$
    – xzczd ♦
    Commented Jan 28, 2023 at 2:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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