# 1) a) Do this for h=0.01 and print the results as requested above.
> 
> restart;
> x0:=0.0; xn:=1.0; y0:=1.0; n:=100; h:=(xn-x0)/n;

                               x0 := 0.


                              xn := 1.0


                              y0 := 1.0


                               n := 100


                          h := .01000000000

> f:=(x,y)->y;    #this is the function f(x,y) in dy/dx = f(x,y)

                           f := (x, y) -> y

> x:=x0; y:=y0;   #initialize x,y for the do loop

                               x := 0.


                               y := 1.0

> for i from 1 to n do
>           k:= f(x,y):  #current slope, use : to suppress output
>           y:= y + h*k: #new y value via Euler
>           x:= x + h:   #updated x-value:
>           if (frac(i/10)=0) then print(x,y,exp(x)) end if
> #         print (x,y,exp(x)):
> od: 

                .1000000000, 1.104622126, 1.105170918


                .2000000000, 1.220190040, 1.221402758


                .3000000000, 1.347848915, 1.349858808


                .4000000000, 1.488863734, 1.491824698


                .5000000000, 1.644631822, 1.648721271


                .6000000000, 1.816696698, 1.822118800


                .7000000000, 2.006763369, 2.013752707


                .8000000000, 2.216715219, 2.225540928


                .9000000000, 2.448632677, 2.459603111


                1.000000000, 2.704813833, 2.718281828

> 
> 
#    b) Print out the graphs for h=0.01 as in figure 2.4.3 of the text.
> restart: with(plots): with(linalg):
Warning, the name changecoords has been redefined

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

> # first we generate the data set:
> x0:=0.0; xn:=1.0; y0:=1.0; n:=100; h:=(xn-x0)/n;

                               x0 := 0.


                              xn := 1.0


                              y0 := 1.0


                               n := 100


                          h := .01000000000

> f:=(x,y)->y;    #this is the function f(x,y) in dy/dx = f(x,y)

                           f := (x, y) -> y

> x:=x0; y:=y0;   #initialize x,y for the do loop

                               x := 0.


                               y := 1.0

> xval:=vector(n+1);yval:=vector(n+1);

                     xval := array(1 .. 101, [])


                     yval := array(1 .. 101, [])

> xval[1]:=x0; yval[1]:=y0; 

                            xval[1] := 0.


                            yval[1] := 1.0

> for i from 1 to n do
>           x:=xval[i]:  #current x
>           y:=yval[i]:  #current y
>           k:= f(x,y):  #current slope
>           yval[i+1]:= y + h*k:   #new y value via Euler
>           xval[i+1]:= x + h:     #updated x-value:
>               od:    
> 
> # approxsol1 is just the points plotted:
> # approxsol2 is the points connected
> approxsol1:=plot([seq([xval[i],yval[i]],
> i=1..n+1)],style=point,color=red, symbol=circle):
> approxsol2:=plot([seq([xval[i],yval[i]], i=1..n+1)],style=line,
> color=red):
> exactsol:=plot(exp(t),t=0..1,color=black):   
> display(approxsol2,exactsol);

> 
> 
> 
> 
> 
#    c) How close is your estimate to e?
> 
> # just look and see
> # notice that if you look closely, you will see that the estimated
> # soution is quite close to the actual solution - BUT is is still
> visibly off!
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
# 2) a) Use improved Euler with n=5 and then n=100 to 
#       see how much better you do for the same initial value problem.  
> 
> 
> f:= (x,y)->y;

                           f := (x, y) -> y

> x0:=0.0; y0:=1.0; xn:=1.0;

                               x0 := 0.


                              y0 := 1.0


                              xn := 1.0

> x:=x0; y:=y0; n:=5.0; h:=(xn-x0)/n;

                               x := 0.


                               y := 1.0


                               n := 5.0


                           h := .2000000000

> for i from 1 to n do
>     k1:=f(x,y):          #left-hand slope
>     k2:=f(x+h,y+h*k1):   #approximation to right-hand slope
>     k:= (k1+k2)/2:       #approximation to average slope
>     y:= y+h*k:            #improved Euler update
>     x:= x+h:             #update x
>     print (x,y,exp(x));
> od:  

                .2000000000, 1.220000000, 1.221402758


                .4000000000, 1.488400000, 1.491824698


                .6000000000, 1.815848000, 1.822118800


                .8000000000, 2.215334560, 2.225540928


                1.000000000, 2.702708163, 2.718281828

# 
> x:=x0; y:=y0; n:=100.0; h:=(xn-x0)/n;

                               x := 0.


                               y := 1.0


                              n := 100.0


                          h := .01000000000

> for i from 1 to n do
>     k1:=f(x,y):          #left-hand slope
>     k2:=f(x+h,y+h*k1):   #approximation to right-hand slope
>     k:= (k1+k2)/2:       #approximation to average slope
>     y:= y+h*k:            #improved Euler update
>     x:= x+h:             #update x
>     if frac(i/10)=0 then print (x,y,exp(x)) end if
> od:

                .1000000000, 1.105169090, 1.105170918


                .2000000000, 1.221398717, 1.221402758


                .3000000000, 1.349852108, 1.349858808


                .4000000000, 1.491814826, 1.491824698


                .5000000000, 1.648707633, 1.648721271


                .6000000000, 1.822100714, 1.822118800


                .7000000000, 2.013729387, 2.013752707


                .8000000000, 2.225511473, 2.225540928


                .9000000000, 2.459566490, 2.459603111


                1.000000000, 2.718236859, 2.718281828

> # much better then before!
> 
# 3) Do Runge-Kutta for the same initial value problem, with n=5 and
#    n=10.  Compare to your results using Euler and improved Euler. 
#     (Again, only do the tables, don't graph.).  
# 
> 
> restart;
> f:= (x,y)->y; x0:=0.0; y0:=1.0; xn:=1.0;

                           f := (x, y) -> y


                               x0 := 0.


                              y0 := 1.0


                              xn := 1.0

> x:=x0; y:=y0; n:=5; h:=(xn-x0)/n;

                               x := 0.


                               y := 1.0


                                n := 5


                           h := .2000000000

> for i from 1 to n do
>     k1:=f(x,y):          #left-hand slope
>     k2:=f(x+h,y+h*k1):   #approximation to right-hand slope
>     k:= (k1+k2)/2:       #approximation to average slope
>     y:= y+h*k:            #improved Euler update
>     x:= x+h:             #update x
>     print(x,y,exp(x));
> od:

                .2000000000, 1.220000000, 1.221402758


                .4000000000, 1.488400000, 1.491824698


                .6000000000, 1.815848000, 1.822118800


                .8000000000, 2.215334560, 2.225540928


                1.000000000, 2.702708163, 2.718281828

> # this is a great approximation with only 10 intervals!
> 
> # here is the n=10 case:
> x:=x0; y:=y0; n:=10; h:=(xn-x0)/n;

                               x := 0.


                               y := 1.0


                               n := 10


                           h := .1000000000

> for i from 1 to n do
>     k1:=f(x,y):          #left-hand slope
>     k2:=f(x+h,y+h*k1):   #approximation to right-hand slope
>     k:= (k1+k2)/2:       #approximation to average slope
>     y:= y+h*k:            #improved Euler update
>     x:= x+h:             #update x
>     print(x,y,exp(x));
> od:

                .1000000000, 1.105000000, 1.105170918


                .2000000000, 1.221025000, 1.221402758


                .3000000000, 1.349232625, 1.349858808


                .4000000000, 1.490902051, 1.491824698


                .5000000000, 1.647446766, 1.648721271


                .6000000000, 1.820428676, 1.822118800


                .7000000000, 2.011573687, 2.013752707


                .8000000000, 2.222788924, 2.225540928


                .9000000000, 2.456181761, 2.459603111


                1.000000000, 2.714080846, 2.718281828

> 
> 
> # a great approximation with only 10 intervals!  Even better then
> n=100 and Euler's method!
> 
# 4) Section 2.4 page 114, #25.
> 
> #a)
> restart:with(DEtools):
> deqn:= 7*x*diff(y(x),x)+y(x)=0;

                               /d      \
                   deqn := 7 x |-- y(x)| + y(x) = 0
                               \dx     /

> dsolve({deqn,y(-1)=1}, y(x));

                                       1/7
                                   (-1)
                            y(x) = -------
                                     1/7
                                    x

> exact:= x-> -1/surd(x,7);

                                          1
                      exact := x -> - ----------
                                      surd(x, 7)

> # b:
> with(plots): with(linalg):
Warning, the name changecoords has been redefined

Warning, the name adjoint has been redefined

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

> # first we generate the data set:
> x0:=-1.0; xn:=0.5; y0:=1.0; h:=.15; n:=ceil((xn-x0)/h);

                              x0 := -1.0


                               xn := .5


                              y0 := 1.0


                               h := .15


                               n := 10

> f:=(x,y)->-y/(7*x);    #this is the function f(x,y) in dy/dx = f(x,y)

                       f := (x, y) -> - 1/7 y/x

> x:=x0; y:=y0;   #initialize x,y for the do loop

                              x := -1.0


                               y := 1.0

> xval:=vector(n+1);   yval:=vector(n+1);

                      xval := array(1 .. 11, [])


                      yval := array(1 .. 11, [])

> xval[1]:=x0; yval[1]:=y0; 

                           xval[1] := -1.0


                            yval[1] := 1.0

> for i from 1 to n do
>           x:=xval[i]:  #current x
>           y:=yval[i]:  #current y
>           k:= f(x,y):  #current slope
>           yval[i+1]:= y + h*k:   #new y value via Euler
>           xval[i+1]:= x + h:     #updated x-value:
>           print(xval[i+1],yval[i+1]);
>               od:    

                          -.85, 1.021428571


                          -.70, 1.047178871


                          -.55, 1.079235367


                          -.40, 1.121283498


                          -.25, 1.181352257


                          -.10, 1.282611022


                           .05, 1.557456241


                           .20, .8899749949


                           .35, .7946205312


                           .50, .7459702946

> 
> # approxsol1 is just the points plotted:
> # approxsol2 is the points connected
> approxsol1:=plot([seq([xval[i],yval[i]],
> i=1..n+1)],style=point,color=red, symbol=circle):
> approxsol2:=plot([seq([xval[i],yval[i]], i=1..n+1)],style=line,
> color=red):
> exactsol:=plot(exact(t),t=-1..1,color=black):
> display(approxsol2,exactsol);

> # it is diffcult to tell from the data alone that there is aproblem.
> # But looking at the graphs and looking at the ODE, it is clear there
> is a 
> # possible problem at x=0
> 
> 
> # c
> # first we do h=0.03
> h:=0.03; n:=ceil((xn-x0)/h);

                               h := .03


                               n := 50

> x:=x0; y:=y0;   #initialize x,y for the do loop

                              x := -1.0


                               y := 1.0

> xval:=vector(n+1);   yval:=vector(n+1);

                      xval := array(1 .. 51, [])


                      yval := array(1 .. 51, [])

> xval[1]:=x0; yval[1]:=y0; 

                           xval[1] := -1.0


                            yval[1] := 1.0

> for i from 1 to n do
>           x:=xval[i]:  #current x
>           y:=yval[i]:  #current y
>           k:= f(x,y):  #current slope
>           yval[i+1]:= y + h*k:   #new y value via Euler
>           xval[i+1]:= x + h:     #updated x-value:
>           if frac((i-1)/5)=0 then print(x,y) end if;
>               od:    

                              -1.0, 1.0


                          -.85, 1.023052521


                          -.70, 1.051187281


                          -.55, 1.087017348


                          -.40, 1.135765593


                          -.25, 1.210371466


                          -.10, 1.361235004


                           .05, 1.872127678


                           .20, 1.471142670


                           .35, 1.350636215

> 
> # approxsol1 is just the points plotted:
> # approxsol2 is the points connected
> approxsol1:=plot([seq([xval[i],yval[i]],
> i=1..n+1)],style=point,color=red, symbol=circle):
> approxsol2:=plot([seq([xval[i],yval[i]], i=1..n+1)],style=line,
> color=red):
> exactsol:=plot(exact(t),t=-1..1,color=black):
> display(approxsol2,exactsol);

> 
> # now we do h=0.006
> h:=0.006; n:=ceil((xn-x0)/h);

                              h := .006


                               n := 250

> x:=x0; y:=y0;   #initialize x,y for the do loop

                              x := -1.0


                               y := 1.0

> xval:=vector(n+1);   yval:=vector(n+1);

                     xval := array(1 .. 251, [])


                     yval := array(1 .. 251, [])

> xval[1]:=x0; yval[1]:=y0; 

                           xval[1] := -1.0


                            yval[1] := 1.0

> for i from 1 to n do
>           x:=xval[i]:  #current x
>           y:=yval[i]:  #current y
>           k:= f(x,y):  #current slope
>           yval[i+1]:= y + h*k:   #new y value via Euler
>           xval[i+1]:= x + h:     #updated x-value:
>           if frac((i-1)/25)=0 then print(x,y) end if;
>               od:    

                              -1.0, 1.0


                          -.850, 1.023400387


                          -.700, 1.052053852


                          -.550, 1.088723590


                          -.400, 1.139018901


                          -.250, 1.217235225


                          -.100, 1.383469420


                          .050, 1.008471774


                          .200, .8210324540


                          .350, .7571469171

> 
> # approxsol1 is just the points plotted:
> # approxsol2 is the points connected
> approxsol1:=plot([seq([xval[i],yval[i]],
> i=1..n+1)],style=point,color=red, symbol=circle):
> approxsol2:=plot([seq([xval[i],yval[i]], i=1..n+1)],style=line,
> color=red):
> exactsol:=plot(exact(t),t=-1..1,color=black):
> display(approxsol2,exactsol);

> # still, the data alone is not convincing (at least to me).
> # but the graph certainly shows that there is some sort of problem
> when x=0
> 
> 
> 
> 
# 5) Section 2.6 page 123:  Work through and understand example #4.
#    Specifically do the following:
#    a) Recreate the data in the table (figure 2.6.9) using your
#       Runge-Kutta routine.  
> 
> restart; with(plots): with(linalg):
Warning, the name changecoords has been redefined

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

> exact:= x->exp(-x);

                        exact := x -> exp(-x)

> f:= (x,y)->5*y-6*exp(-x); x0:=0.0; y0:=1.0; xn:=4.0;

                    f := (x, y) -> 5 y - 6 exp(-x)


                               x0 := 0.


                              y0 := 1.0


                              xn := 4.0

> x:=x0; y:=y0; h:=0.2; n:= ceil((xn-x0)/h);

                               x := 0.


                               y := 1.0


                               h := .2


                               n := 20

> xval:=vector(n+1);   yval:=vector(n+1);

                      xval := array(1 .. 21, [])


                      yval := array(1 .. 21, [])

> xval[1]:=x0; yval[1]:=y0; 

                            xval[1] := 0.


                            yval[1] := 1.0

> for i from 1 to n do
>     x:=xval[i]:
>     y:=yval[i]:
>     k1:=f(x,y):          #left-hand slope
>     k2:=f(x+h/2,y+h*k1/2):   #approximation to right-hand slope
>     k3:=f(x+h/2,y+h*k2/2):
>     k4:=f(x+h, y+h*k3):
> 
>     yval[i+1]:= y+h/6*(k1+2*k2+2*k3+k4):
>     xval[i+1]:= x+h:             #update x
>     if frac(i/2)=0 then print(xval[i+1],yval[i+1]) end if;
> od:

                           .4, .6687956003


                           .8, .4371251820


                           1.2, .2109936451


                          1.6, -.4601900574


                          2.0, -4.721423508


                          2.4, -35.53415136


                          2.8, -261.2502186


                          3.2, -1916.693873


                          3.6, -14059.35439


                          4.0, -103126.5229

> 
> 
> 
> # now for h=0.1
> 
> x:=x0; y:=y0; h:=0.1; n:= ceil((xn-x0)/h);

                               x := 0.


                               y := 1.0


                               h := .1


                               n := 40

> xval:=vector(n+1);   yval:=vector(n+1);

                      xval := array(1 .. 41, [])


                      yval := array(1 .. 41, [])

> xval[1]:=x0; yval[1]:=y0; 

                            xval[1] := 0.


                            yval[1] := 1.0

> for i from 1 to n do
>     x:=xval[i]:
>     y:=yval[i]:
>     k1:=f(x,y):          #left-hand slope
>     k2:=f(x+h/2,y+h*k1/2):   #approximation to right-hand slope
>     k3:=f(x+h/2,y+h*k2/2):
>     k4:=f(x+h, y+h*k3):
> 
>     yval[i+1]:= y+h/6*(k1+2*k2+2*k3+k4):
>     xval[i+1]:= x+h:             #update x
>     if frac(i/4)=0 then print(xval[i+1],yval[i+1]) end if;
> od:

                           .4, .6701959941


                           .8, .4483298137


                           1.2, .2937607749


                           1.6, .1469708761


                          2.0, -.2702590726


                          2.4, -2.904195523


                          2.8, -22.05355342


                          3.2, -163.2510490


                          3.6, -1205.714570


                          4.0, -8903.144048

> 
> 
> 
> 
> # now for h=0.05
> 
> x:=x0; y:=y0; h:=0.05; n:= ceil((xn-x0)/h);

                               x := 0.


                               y := 1.0


                               h := .05


                               n := 80

> xval:=vector(n+1);   yval:=vector(n+1);

                      xval := array(1 .. 81, [])


                      yval := array(1 .. 81, [])

> xval[1]:=x0; yval[1]:=y0; 

                            xval[1] := 0.


                            yval[1] := 1.0

> for i from 1 to n do
>     x:=xval[i]:
>     y:=yval[i]:
>     k1:=f(x,y):          #left-hand slope
>     k2:=f(x+h/2,y+h*k1/2):   #approximation to right-hand slope
>     k3:=f(x+h/2,y+h*k2/2):
>     k4:=f(x+h, y+h*k3):
> 
>     yval[i+1]:= y+h/6*(k1+2*k2+2*k3+k4):
>     xval[i+1]:= x+h:             #update x
>     if frac(i/8)=0 then print(xval[i+1],yval[i+1]) end if;
> od:

                           .40, .6703113029


                           .80, .4492585034


                          1.20, .3006696723


                          1.60, .1980182369


                          2.00, .1066781969


                          2.40, -.1210208487


                          2.80, -1.503657864


                          3.20, -11.51856817


                          3.60, -85.38069959


                          4.00, -631.0329806

> 
> # in this case, the data clearly shows that there is a problem!
> 
> 
>