Euler's method for systems

This worksheet gives code for solving a 2x2 system of ODEs using Maple. Examples illustrate the syntax.

> eulers_method_2x2_exact:=proc(f1,f2,x0,y10,y20,h,n::integer)
local a,i,x,y1,y2,y1n,y2n;
x:=x0;
y1:=y10;
y2:=y20;
a:=[[x,y1,y2]];
for i from 1 to n do
y1n:=y1+h*f1(x,y1,y2);
y2n:=y2+h*f2(x,y1,y2);
x:=x+h;
y1:=y1n;
y2:=y2n;
a:=[op(a),[x,y1,y2]];
od;
RETURN(a);
end;

eulers_method_2x2_exact := proc (f1, f2, x0, y10, y20, h, n::integer) local a, i, x, y1, y2, y1n, y2n; x := x0; y1 := y10; y2 := y20; a := [[x, y1, y2]]; for i to n do y1n := y1+h*f1(x, y1, y2); y2n :...eulers_method_2x2_exact := proc (f1, f2, x0, y10, y20, h, n::integer) local a, i, x, y1, y2, y1n, y2n; x := x0; y1 := y10; y2 := y20; a := [[x, y1, y2]]; for i to n do y1n := y1+h*f1(x, y1, y2); y2n :...eulers_method_2x2_exact := proc (f1, f2, x0, y10, y20, h, n::integer) local a, i, x, y1, y2, y1n, y2n; x := x0; y1 := y10; y2 := y20; a := [[x, y1, y2]]; for i to n do y1n := y1+h*f1(x, y1, y2); y2n :...eulers_method_2x2_exact := proc (f1, f2, x0, y10, y20, h, n::integer) local a, i, x, y1, y2, y1n, y2n; x := x0; y1 := y10; y2 := y20; a := [[x, y1, y2]]; for i to n do y1n := y1+h*f1(x, y1, y2); y2n :...eulers_method_2x2_exact := proc (f1, f2, x0, y10, y20, h, n::integer) local a, i, x, y1, y2, y1n, y2n; x := x0; y1 := y10; y2 := y20; a := [[x, y1, y2]]; for i to n do y1n := y1+h*f1(x, y1, y2); y2n :...eulers_method_2x2_exact := proc (f1, f2, x0, y10, y20, h, n::integer) local a, i, x, y1, y2, y1n, y2n; x := x0; y1 := y10; y2 := y20; a := [[x, y1, y2]]; for i to n do y1n := y1+h*f1(x, y1, y2); y2n :...eulers_method_2x2_exact := proc (f1, f2, x0, y10, y20, h, n::integer) local a, i, x, y1, y2, y1n, y2n; x := x0; y1 := y10; y2 := y20; a := [[x, y1, y2]]; for i to n do y1n := y1+h*f1(x, y1, y2); y2n :...eulers_method_2x2_exact := proc (f1, f2, x0, y10, y20, h, n::integer) local a, i, x, y1, y2, y1n, y2n; x := x0; y1 := y10; y2 := y20; a := [[x, y1, y2]]; for i to n do y1n := y1+h*f1(x, y1, y2); y2n :...eulers_method_2x2_exact := proc (f1, f2, x0, y10, y20, h, n::integer) local a, i, x, y1, y2, y1n, y2n; x := x0; y1 := y10; y2 := y20; a := [[x, y1, y2]]; for i to n do y1n := y1+h*f1(x, y1, y2); y2n :...

> eulers_method_2x2:=proc(f1,f2,x0,y10,y20,h,n::integer)
local a,i,x,y1,y2,y1n,y2n;
x:=evalf(x0);
y1:=evalf(y10);
y2:=evalf(y20);
a:=[[x,y1,y2]];
for i from 1 to n do
y1n:=y1+h*f1(x,y1,y2);
y2n:=y2+h*f2(x,y1,y2);
x:=x+h;
y1:=y1n;
y2:=y2n;
a:=[op(a),[x,y1,y2]];
od;
RETURN(a);
end;

eulers_method_2x2 := proc (f1, f2, x0, y10, y20, h, n::integer) local a, i, x, y1, y2, y1n, y2n; x := evalf(x0); y1 := evalf(y10); y2 := evalf(y20); a := [[x, y1, y2]]; for i to n do y1n := y1+h*f1(x,...eulers_method_2x2 := proc (f1, f2, x0, y10, y20, h, n::integer) local a, i, x, y1, y2, y1n, y2n; x := evalf(x0); y1 := evalf(y10); y2 := evalf(y20); a := [[x, y1, y2]]; for i to n do y1n := y1+h*f1(x,...eulers_method_2x2 := proc (f1, f2, x0, y10, y20, h, n::integer) local a, i, x, y1, y2, y1n, y2n; x := evalf(x0); y1 := evalf(y10); y2 := evalf(y20); a := [[x, y1, y2]]; for i to n do y1n := y1+h*f1(x,...eulers_method_2x2 := proc (f1, f2, x0, y10, y20, h, n::integer) local a, i, x, y1, y2, y1n, y2n; x := evalf(x0); y1 := evalf(y10); y2 := evalf(y20); a := [[x, y1, y2]]; for i to n do y1n := y1+h*f1(x,...eulers_method_2x2 := proc (f1, f2, x0, y10, y20, h, n::integer) local a, i, x, y1, y2, y1n, y2n; x := evalf(x0); y1 := evalf(y10); y2 := evalf(y20); a := [[x, y1, y2]]; for i to n do y1n := y1+h*f1(x,...eulers_method_2x2 := proc (f1, f2, x0, y10, y20, h, n::integer) local a, i, x, y1, y2, y1n, y2n; x := evalf(x0); y1 := evalf(y10); y2 := evalf(y20); a := [[x, y1, y2]]; for i to n do y1n := y1+h*f1(x,...eulers_method_2x2 := proc (f1, f2, x0, y10, y20, h, n::integer) local a, i, x, y1, y2, y1n, y2n; x := evalf(x0); y1 := evalf(y10); y2 := evalf(y20); a := [[x, y1, y2]]; for i to n do y1n := y1+h*f1(x,...eulers_method_2x2 := proc (f1, f2, x0, y10, y20, h, n::integer) local a, i, x, y1, y2, y1n, y2n; x := evalf(x0); y1 := evalf(y10); y2 := evalf(y20); a := [[x, y1, y2]]; for i to n do y1n := y1+h*f1(x,...eulers_method_2x2 := proc (f1, f2, x0, y10, y20, h, n::integer) local a, i, x, y1, y2, y1n, y2n; x := evalf(x0); y1 := evalf(y10); y2 := evalf(y20); a := [[x, y1, y2]]; for i to n do y1n := y1+h*f1(x,...

>

The example below shows how to use n=3 steps of Eulers method with step-size h=1/3

to estimate y1(1) amd y2(1), where y1 and y2 satisfy the system of ODEs

 y1' = y1 - y2 + 1,    y1(0) = 0,   

 y2' = y1 + 3y2 - 1,  y2(0)=0.

So, x0 = 0, y10 = 0, y20 = 0.

> f1:=(x,y1,y2)->y1-y2+1;
f2:=(x,y1,y2)->y1+3*y2-1;
eulers_method_2x2_exact(f1,f2,0,0,0,1/3,3);
eulers_method_2x2(f1,f2,0,0,0,1/3,3);

f1 := proc (x, y1, y2) options operator, arrow; y1-y2+1 end proc

f2 := proc (x, y1, y2) options operator, arrow; y1+3*y2-1 end proc

[[0, 0, 0], [1/3, 1/3, (-1)/3], [2/3, 8/9, (-8)/9], [1, 49/27, (-49)/27]]

[[0., 0., 0.], [.3333333333, .3333333333, -.3333333333], [.6666666666, .8888888890, -.8888888890], [.9999999999, 1.814814815, -1.814814815]]

>

>

>

> sys4 := [diff(x1(t),t) = x2(t), diff(x2(t),t) = 1-2*x1(t)+3*x2(t)];
sol4 := dsolve(sys4);

sys4 := [diff(x1(t), t) = x2(t), diff(x2(t), t) = 1-2*x1(t)+3*x2(t)]

sol4 := {x2(t) = 2*exp(2*t)*_C1+exp(t)*_C2, x1(t) = 1/2+exp(2*t)*_C1+exp(t)*_C2}

> dsolve({diff(x(t),t,t)-3*diff(x(t),t)+2*x(t),x(0)=0,D(x)(0)=1},x(t));

x(t) = -exp(t)+exp(2*t)

>

The example below shows how to use n=4 steps of Eulers method with step-size h=1/4
to estimate y1(1) amd y2(1), where y1 and y2 satisfy the system of ODEs
 y1' = y2,    y1(0) = 0,   
 y2' = -2y1 + 3y2 + 1,  y2(0)=1.
So, x0 = 0, y10 = 0, y20 = 1.

> f1:=(x,y1,y2)->y2;
f2:=(x,y1,y2)->1-2*y1+3*y2;
eulers_method_2x2_exact(f1,f2,0,0,1,1/4,4);
eulers_method_2x2(f1,f2,0,0,1,1/4,4);

f1 := proc (x, y1, y2) options operator, arrow; y2 end proc

f2 := proc (x, y1, y2) options operator, arrow; 1-2*y1+3*y2 end proc

[[0, 0, 1], [1/4, 1/4, 2], [1/2, 3/4, 29/8], [3/4, 53/32, 199/32], [1, 411/128, 1319/128]]

[[0., 0., 1.], [.2500000000, .2500000000, 2.000000000], [.5000000000, .7500000000, 3.625000000], [.7500000000, 1.656250000, 6.218750000], [1.000000000, 3.210937500, 10.30468750]][[0., 0., 1.], [.2500000000, .2500000000, 2.000000000], [.5000000000, .7500000000, 3.625000000], [.7500000000, 1.656250000, 6.218750000], [1.000000000, 3.210937500, 10.30468750]]

>

Last modified 7-28-2003 by wdj.