# SECTION 5.3

# Example 2

with(linalg);

Warning, the protected names norm and trace have been redefined and unprotected

NiM3XnIlLkJsb2NrRGlhZ29uYWxHJSxHcmFtU2NobWlkdEclLEpvcmRhbkJsb2NrRyUpTFVkZWNvbXBHJSlRUmRlY29tcEclKldyb25za2lhbkclJ2FkZGNvbEclJ2FkZHJvd0clJGFkakclKGFkam9pbnRHJSZhbmdsZUclKGF1Z21lbnRHJShiYWNrc3ViRyUlYmFuZEclJmJhc2lzRyUnYmV6b3V0RyUsYmxvY2ttYXRyaXhHJShjaGFybWF0RyUpY2hhcnBvbHlHJSljaG9sZXNreUclJGNvbEclJ2NvbGRpbUclKWNvbHNwYWNlRyUoY29sc3BhbkclKmNvbXBhbmlvbkclJ2NvbmNhdEclJWNvbmRHJSljb3B5aW50b0clKmNyb3NzcHJvZEclJWN1cmxHJSlkZWZpbml0ZUclKGRlbGNvbHNHJShkZWxyb3dzRyUkZGV0RyUlZGlhZ0clKGRpdmVyZ2VHJShkb3Rwcm9kRyUqZWlnZW52YWxzRyUsZWlnZW52YWx1ZXNHJS1laWdlbnZlY3RvcnNHJStlaWdlbnZlY3RzRyUsZW50ZXJtYXRyaXhHJSZlcXVhbEclLGV4cG9uZW50aWFsRyUnZXh0ZW5kRyUsZmZnYXVzc2VsaW1HJSpmaWJvbmFjY2lHJStmb3J3YXJkc3ViRyUqZnJvYmVuaXVzRyUqZ2F1c3NlbGltRyUqZ2F1c3Nqb3JkRyUoZ2VuZXFuc0clKmdlbm1hdHJpeEclJWdyYWRHJSloYWRhbWFyZEclKGhlcm1pdGVHJShoZXNzaWFuRyUoaGlsYmVydEclK2h0cmFuc3Bvc2VHJSlpaGVybWl0ZUclKmluZGV4ZnVuY0clKmlubmVycHJvZEclKWludGJhc2lzRyUoaW52ZXJzZUclJ2lzbWl0aEclKmlzc2ltaWxhckclJ2lzemVyb0clKWphY29iaWFuRyUnam9yZGFuRyUna2VybmVsRyUqbGFwbGFjaWFuRyUqbGVhc3RzcXJzRyUpbGluc29sdmVHJSdtYXRhZGRHJSdtYXRyaXhHJSZtaW5vckclKG1pbnBvbHlHJSdtdWxjb2xHJSdtdWxyb3dHJSltdWx0aXBseUclJW5vcm1HJSpub3JtYWxpemVHJSpudWxsc3BhY2VHJSdvcnRob2dHJSpwZXJtYW5lbnRHJSZwaXZvdEclKnBvdGVudGlhbEclK3JhbmRtYXRyaXhHJStyYW5kdmVjdG9yRyUlcmFua0clKHJhdGZvcm1HJSRyb3dHJSdyb3dkaW1HJSlyb3dzcGFjZUclKHJvd3NwYW5HJSVycmVmRyUqc2NhbGFybXVsRyUtc2luZ3VsYXJ2YWxzRyUmc21pdGhHJSxzdGFja21hdHJpeEclKnN1Ym1hdHJpeEclKnN1YnZlY3RvckclKXN1bWJhc2lzRyUoc3dhcGNvbEclKHN3YXByb3dHJSpzeWx2ZXN0ZXJHJSl0b2VwbGl0ekclJnRyYWNlRyUqdHJhbnNwb3NlRyUsdmFuZGVybW9uZGVHJSp2ZWNwb3RlbnRHJSh2ZWN0ZGltRyUndmVjdG9yRyUqd3JvbnNraWFuRw==

x := vector(2);

NiM+JSJ4Ry0lJmFycmF5RzYkOyIiIiIiIzci

f := vector(2);

NiM+JSJmRy0lJmFycmF5RzYkOyIiIiIiIzci

f := x -> [x[2],-4*x[2]-3*x[1]];

NiM+JSJmR2orNiMlInhHNiI2JCUpb3BlcmF0b3JHJSZhcnJvd0dGKDckJjkkNiMiIiMsJiomIiIlIiIiRi1GNCEiIiomIiIkRjQmRi42I0Y0RjRGNUYoRihGKDYjIiIh

x0 := vector(2); xnew := vector(2);

NiM+JSN4MEctJSZhcnJheUc2JDsiIiIiIiM3Ig==

NiM+JSV4bmV3Ry0lJmFycmF5RzYkOyIiIiIiIzci

x0[1] := 1.5; x0[2] := -2.5;

NiM+JiUjeDBHNiMiIiIkIiM6ISIi

NiM+JiUjeDBHNiMiIiMkISNEISIi

h := .1;

NiM+JSJoRyQiIiIhIiI=

x := evalm(x0);

NiM+JSJ4R0slJ3ZlY3Rvckc2IzckJCIjOiEiIiQhI0RGK1EocHByaW50MDYi

# First step:

x := evalm(x+h*f(x));

NiM+JSJ4R0slJ3ZlY3Rvckc2IzckJCIkRCIhIiMkISQmPkYrUShwcHJpbnQxNiI=

# Second step:

x := evalm(x + h*f(x));

NiM+JSJ4R0slJ3ZlY3Rvckc2IzckJCIlYjUhIiQkISVYOkYrUShwcHJpbnQyNiI=

# Now take all ten steps:

# First check syntax of "for" loop:

?for

t := 0;x := evalm(x0);# Initialization of t and x

NiM+JSJ0RyIiIQ==

NiM+JSJ4R0slJ3ZlY3Rvckc2IzckJCIjOiEiIiQhI0RGK1EocHByaW50MzYi

for n from 1 to 10 do

t := t + h;x := evalm(x+h*f(x));# Next t and x

end do;

NiM+JSJ0RyQiIiIhIiI=

NiM+JSJ4R0slJ3ZlY3Rvckc2IzckJCIkRCIhIiMkISQmPkYrUShwcHJpbnQ0NiI=

NiM+JSJ0RyQiIiMhIiI=

NiM+JSJ4R0slJ3ZlY3Rvckc2IzckJCIlYjUhIiQkISVYOkYrUShwcHJpbnQ1NiI=

NiM+JSJ0RyQiIiQhIiI=

NiM+JSJ4R0slJ3ZlY3Rvckc2IzckJCIlMCEqISIlJCEmTkMiRitRKHBwcmludDY2Ig==

NiM+JSJ0RyQiIiUhIiI=

NiM+JSJ4R0slJ3ZlY3Rvckc2IzckJCImOncoISImJCEnRDs1RitRKHBwcmludDc2Ig==

NiM+JSJ0RyQiIiYhIiI=

NiM+JSJ4R0slJ3ZlY3Rvckc2IzckJCInRFhuISInJCEnJmZVKUYrUShwcHJpbnQ4NiI=

NiM+JSJ0RyQiIichIiI=

NiM+JSJ4R0slJ3ZlY3Rvckc2IzckJCIoYkUhZiEiKCQhKFgienFGK1EocHByaW50OTYi

NiM+JSJ0RyQiIighIiI=

NiM+JSJ4R0slJ3ZlY3Rvckc2IzckJCIpMHUlPiYhIikkISlORz1nRitRKXBwcmludDEwNiI=

NiM+JSJ0RyQiIikhIiI=

NiM+JSJ4R0slJ3ZlY3Rvckc2IzckJCIqOjdIZiUhIiokISpEI1JwXkYrUSlwcHJpbnQxMTYi

NiM+JSJ0RyQiIiohIiI=

NiM+JSJ4R0slJ3ZlY3Rvckc2IzckJCIrREgoZjIlISM1JCErJioqMyZ6V0YrUSlwcHJpbnQxMjYi

NiM+JSJ0RyQiIzUhIiI=

NiM+JSJ4R0slJ3ZlY3Rvckc2IzckJCIrRT8tR08hIzUkISt2c1w1UkYrUSlwcHJpbnQxMzYi

# Example 3

k := 'k'; n := 'n';

NiM+JSJrR0Yk

NiM+JSJuR0Yk

# We will put the above loop inside another loop that changes the step sizes h.

# Also, we'll only display h and the approximation for y at t=1:

print('h','y(1,h)');

NiQlImhHLSUieUc2JCIiIkYj

h := .1: N := 10: # Initial step size and number of steps

for k from 1 to 5 do

t := 0:x := evalm(x0): # Initialization of t and x

for n from 1 to N do

t := t + h: x := evalm(x+h*f(x)):# Next t and x

end do:

print(h,x[1]);

h := h/2: N := 2*N:

end do:

NiQkIiIiISIiJCIrRT8tR08hIzU=

NiQkIisrKysrXSEjNiQiKyIpb2x5UCEjNQ==

NiQkIisrKysrRCEjNiQiK1QsWGBRISM1

NiQkIisrKytdNyEjNiQiK1AvbSEqUSEjNQ==

NiQkIisrKytdaSEjNyQiK0tNQDRSISM1

# Example 4

f := x -> [2*x[1]-2*x[1]*x[2],x[1]*x[2]-x[2]];

NiM+JSJmR2orNiMlInhHNiI2JCUpb3BlcmF0b3JHJSZhcnJvd0dGKDckLCYqJiIiIyIiIiY5JDYjRjBGMEYwKihGL0YwRjFGMCZGMjYjRi9GMCEiIiwmKiZGMUYwRjVGMEYwRjVGN0YoRihGKDYjIiIh

x0 := [1,3];

NiM+JSN4MEc3JCIiIiIiJA==

h := 1.; x := evalm(x0); # Initialization

NiM+JSJoRyQiIiIiIiE=

NiM+JSJ4R0slJ3ZlY3Rvckc2IzckIiIiIiIkUSlwcHJpbnQxNDYi

k1 := evalm(h*f(x));

NiM+JSNrMUdLJSd2ZWN0b3JHNiM3JCQhIiUiIiEkRitGK1EpcHByaW50MTU2Ig==

k2 := evalm(h*f(evalm(x+.5*k1)));

NiM+JSNrMkdLJSd2ZWN0b3JHNiM3JCQiI1MhIiIkISNnRitRKXBwcmludDE2NiI=

k3 := evalm(h*f(evalm(x+.5*k2)));

NiM+JSNrM0dLJSd2ZWN0b3JHNiM3JCQiJCsnISIjJCIiIUYtUSlwcHJpbnQxNzYi

k4 := evalm(h*f(evalm(x+k3)));

NiM+JSNrNEdLJSd2ZWN0b3JHNiM3JCQhJStHISIjJCIlKz1GK1EpcHByaW50MTg2Ig==

x := evalm(x + (k1+2*(k2+k3)+k4)/6.);

NiM+JSJ4R0slJ3ZlY3Rvckc2IzckJCErLCsrKzUhIiokIissKysrU0YrUSlwcHJpbnQxOTYi

# The above represents a single Runge-Kutta step. It will be

# convenient to make it into a "procedure" called rk4. First, we review the appropriate syntax:

? proc

rk4 := proc(x,h)

local k1, k2, k3, k4;

k1 := evalm(h*f(x));

k2 := evalm(h*f(evalm(x+.5*k1)));

k3 := evalm(h*f(evalm(x+.5*k2)));

k4 := evalm(h*f(evalm(x+k3)));

x := evalm(x + (k1+2*(k2+k3)+k4)/6.);

end proc;

NiM+JSRyazRHais2JCUieEclImhHNiYlI2sxRyUjazJHJSNrM0clI2s0RzYiRi5DJz44JC0lJmV2YWxtRzYjKiY5JSIiIi0lImZHNiM5JEY3PjglLUYzNiMqJkY2RjctRjk2Iy1GMzYjLCZGO0Y3KiYkIiImISIiRjdGMUY3RjdGNz44Ji1GMzYjKiZGNkY3LUY5NiMtRjM2IywmRjtGNyomRkdGN0Y9RjdGN0Y3PjgnLUYzNiMqJkY2RjctRjk2Iy1GMzYjLCZGO0Y3RktGN0Y3PkY7LUYzNiMsJkY7RjcqJiwqRjFGNyomIiIjRjdGPUY3RjcqJkZgb0Y3RktGN0Y3RlZGN0Y3JCIiJyIiIUZJRjdGLkYuRi42JEZkb0Zkbw==

# Next, we'll test the procedure to see if it gives the same result that we got before:

x := evalm(x0); h := 1.;

NiM+JSJ4R0slJ3ZlY3Rvckc2IzckIiIiIiIkUSlwcHJpbnQyMDYi

NiM+JSJoRyQiIiIiIiE=

xnew := rk4(x,h);

NiM+JSV4bmV3R0slJ3ZlY3Rvckc2IzckJCErLCsrKzUhIiokIissKysrU0YrUSlwcHJpbnQyMTYi

# Next, we create the loops to change h, and to approximate x at t=1 for each h:

n := 'n'; k := 'k';

NiM+JSJuR0Yk

NiM+JSJrR0Yk

print('h','x(1,h)');

NiQlImhHLSUieEc2JCIiIkYj

h := 1.: N := 1: # Initial step size and number of steps

for k from 1 to 5 do

t := 0:x := evalm(x0): # Initialization of t and x

for n from 1 to N do

x := rk4(x,h):end do:

print(h,x);

h := h/2: N := 2*N:

end do:

NiQkIiIiIiIhSyUndmVjdG9yRzYjNyQkISssKysrNSEiKiQiKywrKytTRixRKXBwcmludDIyNiI=

NiQkIisrKysrXSEjNUslJ3ZlY3Rvckc2IzckJCIrRCcpPW05RiUkIitEL2N0OSEiKlEpcHByaW50MjM2Ig==

NiQkIisrKysrRCEjNUslJ3ZlY3Rvckc2IzckJCIreUVWJil5ISM2JCIrRDlwazkhIipRKXBwcmludDI0NiI=

NiQkIisrKytdNyEjNUslJ3ZlY3Rvckc2IzckJCIrWy1oU3ghIzYkIis6ZVhrOSEiKlEpcHByaW50MjU2Ig==

NiQkIisrKytdaSEjNkslJ3ZlY3Rvckc2IzckJCIrIypmc014RiUkIitEJltXWSIhIipRKXBwcmludDI2NiI=