4

How can I implement this pseudocode in Mathematica?

Algorithm:

(* input: f, x0, xmax, y0, n  *)

h = (xmax-x0)/n; u[1] = x0; v[1] = y0; For i = 2,n+1: u[i] = u[i-1] + h; v[i] = v[i-1] + h * f( u[i-1] , v[i-1] ); End for. Return {{u[1],v[1]},{u[2],v[2]},...,{u[n+1],v[n+1]}};

This is my take on it but it's not working properly.

MyEuler[f, x0, xmax, y0, n] = Module[{},
  h = (xmax - x0)/n;
  u[1] = x0;
  v[1] = y0;
  For[i = 2, i <= n + 2, i++,
    u[i] = u[i - 1] + h;
    v[i] = v[i - 1] + h*f (u[i - 1], v[i - 1]);]
   Table[{u[k], v[k]}, {k, 1, n + 1}]]
Michael E2
  • 235,386
  • 17
  • 334
  • 747
Birgitt
  • 267
  • 1
  • 4
  • 3
    f (u[i - 1], v[i - 1]) should probably be f[u[i - 1], v[i - 1]] if f is a function. And why end your loop at at n+2 instead of n+1 is a difference, but perhaps an unimportant one. Other than that, read a good introductory text on how to program Mathematica. You aren't taking good advantage of its features, other than its quasi-ability to pretend to be C or Java. – Michael E2 Nov 03 '21 at 12:18
  • 3
    Possible duplicates: (22042), (22413*), (41783), (43623), (206106), especially (22413). – Michael E2 Nov 03 '21 at 12:29
  • 2
    Further, arguments of a function are defined by underscore like: `MyEuler[f_,x0_, ....] – Daniel Huber Nov 03 '21 at 13:06
  • 1
    Are you aware of the fact that Euler method is already implemented into the standard NDSolve function and can be fixed by the option Method? It is just to make sure that we understood you correctly. – Alexei Boulbitch Nov 04 '21 at 15:40
  • Yes I am aware,but my teacher specifically told us to solve the problem this way. – Birgitt Nov 04 '21 at 22:18

2 Answers2

6

In a functional style, I recommend NestList, the key trick of which is to properly define an iteration function to be used repeatedly via pure/anonymous Function (&). And more can be found out in @MichaelE2 's comments.

Clear[myEuler]
myEuler[f : (_Symbol | _Function), x0_Real, xmax_Real, y0_Real, n_Integer] :=
Module[{h = (xmax - x0)/n, iterate},
    iterate = # + h {1, f @@ #} &;
    NestList[iterate, {x0, y0}, n]
]

Note that, considering the numeric nature of the problem, here I set explicitly type checking for each argument of myEuler. The consequence is that unless all the type checking returns True, myEuler will refuse to run. In this way would it be disentangled with MMA's symbolic ability to some extend, which is actually not needed here.

As an example, I set $ f(x, y) = x^2 + y^2, x_0 = 0., x_\text{max} = 1., y_0 = .1 $, and $ n = 100 $, and get the result plotted as below:

enter image description here

5

There are several syntactical mistakes in your code.

  1. Function arguments are missing an underscore:

MyEuler[f_, x0_, xmax_, y0_, n_]

  1. You should use SetDelayed instead of Set, that is: you are missing a colon:

MyEuler[f_, x0_, xmax_, y0_, n_] :=

  1. Calling a function requires square brackets, not round ones:

f[u[i - 1], v[i - 1]]

  1. Separate For from Table by the additional semicolon:

For[ ... ]; Table[ ... ]

Put all of this together, and the code should work:

MyEuler[f_, x0_, xmax_, y0_, n_] := Module[{},
  h = (xmax - x0)/n;
  u[1] = x0;
  v[1] = y0;
  For[i = 2, i <= n + 2, i++,
   u[i] = u[i - 1] + h;
   v[i] = v[i - 1] + h*f[u[i - 1], v[i - 1]]
   ];
  Table[{u[k], v[k]}, {k, 1, n + 1}]
  ]

f[x_, y_] := x + y;

MyEuler[f, 0, 1, 0, 3]

(* {{0, 0}, {1/3, 0}, {2/3, 1/9}, {1, 10/27}} *)

Domen
  • 23,608
  • 1
  • 27
  • 45
  • We could put h = (xmax - x0)/n, u, v} inside the currently empty 1st argument of Module. And given that u and v are given initial values named x0 and y0, it would be sensible to change indexing to begin with 0 instead of 1 and to then run i only up to n. – murray Nov 03 '21 at 14:53
  • @murray, yes, certainly there can be several improvements of the code (as pointed out in other answers). But since it looks like the OP is a beginner in Mathematica, I tried to make only the most necessary changes to make the code work. – Domen Nov 03 '21 at 15:10