Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

share|improve this question
    
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 8 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 8 at 12:42

1 Answer 1

up vote 4 down vote accepted

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

share|improve this answer
    
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 8 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 en.wikipedia.org/wiki/Euler_method (particularly the 1 d example $f(t,y(t))=y$ –  ubpdqn Mar 8 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 8 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 8 at 12:58
    
@MichaelE2 sorry for sensitivity...always keen to correct my own errors and misconceptions. –  ubpdqn Mar 8 at 13:03

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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