# SECTION 3.6# Example 1# Remark: This solution is a modification of the solution to Example 2, Sect 1.4f := (x,y) -> y;

NiM+JSJmR2orNiQlInhHJSJ5RzYiNiQlKW9wZXJhdG9yRyUmYXJyb3dHRik5JUYpRilGKTYkIiIhRi8=

x0 := 0; y0 := 1.; e := exp(1.);

NiM+JSN4MEciIiE=NiM+JSN5MEckIiIiIiIhNiM+JSJlRyQiK0c9Rz1GISIq

Numbersteps := [1,10,10^2,10^3,10^4];

NiM+JSxOdW1iZXJzdGVwc0c3JyIiIiIjNSIkKyIiJSs1IiYrKyI=

print('h','Approx','Error','Error/h^2');

NiYlImhHJSdBcHByb3hHJSZFcnJvckcqJkYlIiIiRiMhIiM=

for k from 1 to 5 doN := Numbersteps[k]:h := evalf(1/N):x := x0: y := y0:for n from 1 to N dox1 := x+h: f1 := f(x,y): # Use new variables to avoid repeating calculations y := y + h/2*(f1 + f(x1,y+h*f1)):x := x1:end do:Error := e-y:print(h,y,Error,Error/h^2);end do:

NiYkIiIiIiIhJCIrKysrK0QhIiokIipHPUc9I0YoJCIrIUc9Rz0jISM1NiYkIisrKysrNSEjNSQiK1kzMzlGISIqJCIoIyk0PyVGKCQiKys/KTQ/JUYlNiYkIisrKysrNSEjNiQiK2ZvQj1GISIqJCImcFwlRigkIisrKyFwXCUhIzU=NiYkIisrKysrNSEjNyQiK3o4Rz1GISIqJCIkXCVGKCQiKysrKyFcJSEjNQ==NiYkIisrKysrNSEjOCQiK208Rz1GISIqJCIjaUYoJCIrKysrK2lGKA==

# Remark: The last two results for Error/h^2 differ from the text result because of rounding errors. One can see this by increasing the number of digits used in calculations, as shown below, and then repeating the above example. Then reset Digits to 10.Digits;

NiMiIzU=

Digits := 20;

NiM+JSdEaWdpdHNHIiM/

Digits := 10;

NiM+JSdEaWdpdHNHIiM1

# Example 2# We need to write a subroutine as in the text on page 129, but will not fullyimplement the algorithm on page 130, for simplicity. First, we review the syntax for the "proc" command:?procieuler := proc(x0,y0,c,N,PRNTR)local h,x,y,i,F,G: h := evalf( (c-x0)/N); # We use evalf to make sure h is a decimal number.x := x0: y := y0:for i from 1 to N doF := f(x,y): G := f(x+h,y+h*F):x := x+h:y := y + h*(F+G)/2:if(PRNTR=1) thenprint(x,y);end if:end do:return(x,y):end proc;

NiM+JSdpZXVsZXJHais2JyUjeDBHJSN5MEclImNHJSJORyUmUFJOVFJHNiglImhHJSJ4RyUieUclImlHJSJGRyUiR0c2IkYzQyc+OCQtJSZldmFsZkc2IyomLCY5JiIiIjkkISIiRj05J0Y/PjglRj4+OCY5JT8oOCdGPUY9RkAlJXRydWVHQyc+OCgtJSJmRzYkRkJGRD44KS1GTTYkLCZGQkY9RjZGPSwmRkRGPSomRjZGPUZLRj1GPT5GQkZTPkZELCZGREY9KiYjRj0iIiNGPSomRjZGPSwmRktGPUZQRj1GPUY9Rj1AJC85KEY9LSUmcHJpbnRHRk5PRk5GM0YzRjM2JyIiIUZfb0ZfbyIjdCIqNy5AUyI=

# We test the procedure by using it on the last example, with h =1 and h=.1:# Note that x0, y0, and f have already been defined.ieuler(x0,y0,1,1,0);

NiQkIiIiIiIhJCIrKysrK0QhIio=

ieuler(x0,y0,1,10,0);

NiQkIisrKysrNSEiKiQiK1kzMzlGRiU=

# If we want to save the output from ieuler, we can do the following: out := ieuler(x0,y0,1,10,0);

NiM+JSRvdXRHNiQkIisrKysrNSEiKiQiK1kzMzlGRig=

out[1];out[2];

NiMkIisrKysrNSEiKg==NiMkIitZMzM5RiEiKg==

# Next, we apply ieuler to Example 2:f := (x,y) -> x+2*y;

NiM+JSJmR2orNiQlInhHJSJ5RzYiNiQlKW9wZXJhdG9yRyUmYXJyb3dHRiksJjkkIiIiKiYiIiNGLzklRi9GL0YpRilGKTYkIiIhRjQ=

x0 := 0; y0 := .25; c := 2.;

NiM+JSN4MEciIiE=NiM+JSN5MEckIiNEISIjNiM+JSJjRyQiIiMiIiE=

ieuler(x0,y0,c,1,0);# This is the h=2 case.

NiQkIiIjIiIhJCIrKysrXV8hIio=

# We next let h take the values 2/2^k, k = 1,...,10# Since h = 2/N, this means N =2/h will have the values = 2^kprint('k','h=2/2^k','y(2,h)');

NiUlImtHLyUiaEcsJComIiIjIiIiKUYoRiMhIiJGKS0lInlHNiRGKEYl

for k from 1 to 10 doout := ieuler(x0,y0,c,2^k,0):print(k,2/2^k,out[2]);end do:

NiUiIiJGIyQiKysrK0Q2ISIpNiUiIiMjIiIiRiMkIisrXTdHPSEiKQ==NiUiIiQjIiIiIiIlJCIrIypwMTFCISIpNiUiIiUjIiIiIiIpJCIrNjwsN0QhIik=NiUiIiYjIiIiIiM7JCIrcW43ekQhIik=NiUiIicjIiIiIiNLJCIrNzw4KWYjISIpNiUiIigjIiIiIiNrJCIrLD48LkUhIik=NiUiIikjIiIiIiRHIiQiKyZSb1dnIyEiKQ==NiUiIiojIiIiIiRjIyQiK2Jxei9FISIpNiUiIzUjIiIiIiQ3JiQiKzApeltnIyEiKQ==