2

I need to program Euler's method to solve a system of two diffferential equations of first order.

Fist, I have programmed the Euler's method for just one differential equation:

 euler[f_, ini_, int_, h_] := 
 Module[{x, y, l}, l = {ini}; 
 y = ini[[2]]; 
 For[x = int[[1]], x <= int[[2]] - h, x = x + h, 
 y = y + h*f[x, y]; l = AppendTo[l, {x + h, y}]];l]

And it works fine.

f[x_, y_] := 3 (1 - y);
euler[f, {0, 0.01}, {0, 10}, 1]
{{0, 0.01}, {1, 2.98}, {2, -2.96}, {3,8.92}, {4, -14.84}, {5, 32.68}, {6, -62.36}, {7,127.72}, {8, -252.44}, {9, 507.88}, {10, -1012.76}}

Now, changing a little the previous function I want to solve a system of two differential equations of first order. I programmed the following:

euler2[f_, g_, inif_,inig_, int_, h_] :=  
Module[{x, y, t, list1, list2}, list1 = {inif}; list2 = {inig}; x = inif[[2]]; 
y = inig[[2]]; For[t = int[[1]], t <= int[[2]] - h, t = t + h, 
x = x + h*f[x, y, t]; y = y + h*g[x, y, t]]; 
list1 = AppendTo[list1, {t + h, x}]; 
list2 = AppendTo[list2, {t + h, y}]; {list1, list2}]

But I got an error. What am I doing wrong?

For example, I define f anf g in this way:

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

inif= is a set of the initial value of x (x(0)=1$\to$ inif={0,1}) and inig= is a set of the initial value of y (y(0)=1 $\to$ inig={0,1})

int is the interval where I want to calculate the solution int={0,10} and h the lenght of each step h=1

Thank you for your help.

J. M.'s missing motivation
  • 124,525
  • 11
  • 401
  • 574
Triana
  • 49
  • 1
  • 9
  • f[x_, y_, t_] := y; why would you pass in x and t to the function f if they are not used in the body of the function? g[x_, y_, t_] :=0.5x+0.5y and here, t is not used in the body of the function, so why is it passed? if you have the problem description you are trying to solve from the source, can you post it? I am having hard time seeing what you are trying to do. – Nasser Mar 08 '14 at 11:58
  • @Nasser: It was just an example, but I want to programm the method for any system of two differential equations. I used t because it is the variable of x and y- (x[t], y[t]). I have edited a little. I hope you can understand a little more. – Triana Mar 08 '14 at 12:42

1 Answers1

4

I believe your code is not working as you are not incrementing your iterator t. I present the following (using FoldList) as alternative iterative approach. 

euler2d[f_, g_, h_, t0_, f0_, g0_, n_] :=
 FoldList[#1 + 
    h {f[#2, First@#1, Last@#1], g[#2, First@#1, Last@#1]} &, {f0, 
   g0}, Range[t0, t0 + n h , h]]

Note in this case the interval of interest: [t0,t0+n h]. You could recode as desired.

Test functions:

k[t_, x_, y_] := -y
v[t_, x_, y_] := (x - 2 y)/2

Exact solution:

sol = {x[t], y[t]} /. 
  First@DSolve[{x'[t] == -y[t], y'[t] == (x[t] - 2 y[t])/2, x[0] == 1,
      y[0] == 1}, {x[t], y[t]}, t]

Comparing for interval [0,3]:

Show[ParametricPlot[sol /. t -> u, {u, 0, 3}], 
 ListPlot[euler2d[k, v, 0.05, 0, 1, 1, 60], 
  PlotStyle -> Red], Frame -> True, 
 FrameLabel -> {TraditionalForm@x, TraditionalForm@y}, 
 BaseStyle -> 12]

enter image description here

ubpdqn
  • 60,617
  • 3
  • 59
  • 148
  • k[t_, x_, y_] := -y why is this function defined to accept t and x if they do not appear in the body of the function? v[t_, x_, y_] := (x - 2 y)/2, same here. Where is t is the RHS? – Nasser Mar 08 '14 at 12:39
  • The iterated functions have three arguments. These are merely functions that do not explicitly depend on t (i.e they are autonomous). t is usually time. It just means the system is time invariant. t is just use to parametrize. I used linear systems as Euler is first order and will breakdown for non-linear systems. To be explicit, the t is just where the system is at time t. $y'(t)=f(t,y(t))$ See http://en.wikipedia.org/wiki/Euler_method (particularly the 1 d example $f(t,y(t))=y$ – ubpdqn Mar 08 '14 at 12:53
  • @Nasser I believe the OP has set it up that way because it is the generic Euler method scheme for $dx/dt = f(t,x,y), dy/dt = g(t,x,y)$. Functions from a nonautonomous system can be substituted into the scheme without having to alter the program. – Michael E2 Mar 08 '14 at 12:54
  • @MichaelE2 I merely illustrated an autonomous linear system slightly more interesting than the test case in question (which is a straight line). My code generalizes and have tested on linear , eg, x+t. Just wanted to illustrate advantage of FoldList over For. – ubpdqn Mar 08 '14 at 12:58
  • @MichaelE2 sorry for sensitivity...always keen to correct my own errors and misconceptions. – ubpdqn Mar 08 '14 at 13:03
  • @MichaelE2 thanks...I have corrected. When I originally wrote answer I was thinking of a function that you could choose your own symbols...then realized i was overcomplicating: changed the function definition but used old function definition in visualization...thanks again. should be ok now – ubpdqn Mar 08 '14 at 13:11
  • 1
    @Moo it’s been 8 years since I looked at this. I would suggest you could play with ‘skeleton’ of the FoldList approach and adapt and see. Unfortunately, at present I am poor of time and health. The fun is in trying, however. Good luck. – ubpdqn Mar 27 '22 at 07:39