Math 2250

Homework Due Wednesday September 20

Numerical solutions to differential equations

#5, sections 2.4-2.6: I have saved myself time by copying and pasting pieces of code from the Sept 13 class/internet handout, which you can access via the internet. Of course the code needed minor changes to reflect the fact that we are doing a different initial value problem.

> restart:

> with(DEtools):

> dsolve({diff(y(x),x)=y(x)-x-1,y(0)=1},y(x));
#find the solution to our initial value problem

> sol:= x->x+2-exp(x); #our solution

> x0:=0.0; xn:=0.5; y0:=1.0; n:=5; h:=(xn-x0)/n;
#initial values, 5 steps of step size 0.1.

>

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

Here are the results for step-size = 0.1:

> x:=x0;y:=y0;
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:

> print(x,y,sol(x)); #display current values,

> #and compare to exact solution.

> od: #``od'' ends a do loop

So Euler approximation at x=.5 is about 0.89 vs. actual solution of about 0.85.

Now repeat with h=0.05:

> x:=x0;y:=y0;n:=10;h:=0.05;

> x:=x0;y:=y0;
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/2)=0

> then print(x,y,sol(x)); #display current values
fi; #other time

> #and compare to exact solution.

> od: #``od'' ends a do loop

Now for improved Euler, #5 in section 2.5:

> x:=x0;y:=y0;n:=5;h:=.1;

> 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,sol(x));

> od:

So improved Euler did better than unimproved Euler: its value at x=0.5 is 0.853 vs. the actual value of 0.851, and unimproved Euler's calculation of 0.89. Notice we also did much better than unimproved Euler with h=0.05.

Finally, use Runge-Kutta: (This is #5 from 2.6) Notice we only have to do two steps, with h=0.25,

> x:=x0;y:=y0;n:=2;h:=.25;

> for i from 1 to n do

> k1:=f(x,y): #left-hand slope

> k2:=f(x+h/2,y+h*k1/2): #1st guess at midpoint slope

> k3:=f(x+h/2,y+h*k2/2): #second guess at midpoint slope

> k4:=f(x+h,y+h*k3): #guess at right-hand slope

> k:=(k1+2*k2+2*k3+k4)/6: #Simpson's approximation for the integral

> x:=x+h: #x update

> y:=y+h*k: #y update

> print(x,y,sol(x)); #display current values

> od:

Pretty Good: Runge Kutta rounded to five decimal places 0.85130 whereas the exact solution rounds to 0.85128, and error of 0.00002.

#25 section 2.5:

> restart:

> with(DEtools):

> sol:=t->294.0*exp(-t/25.0) - 245;
#from the book; could also have used dsolve.

> f:=(t,v)->-0.04*v -9.8;
t0:=0;v0:=49;
t:=t0;v:=v0;n:=50;h:=10.0/n; #50 subdivisions

> for i from 1 to n do

> k1:=f(t,v): #left-hand slope

> k2:=f(t+h,v+h*k1): #approximation to right-hand slope

> k:= (k1+k2)/2: #approximation to average slope

> v:= v+h*k: #improved Euler update for v

> t:= t+h: #update t
if frac(i/5)=0

> then print(t,v,sol(t),sol(t)-v);
#also print the difference between exact solution
#and approximate solution.
fi;

> od:

> t:=t0;v:=v0;n:=100;h:=10.0/n; #100 subdivisions

> for i from 1 to n do

> k1:=f(t,v): #left-hand slope

> k2:=f(t+h,v+h*k1): #approximation to right-hand slope

> k:= (k1+k2)/2: #approximation to average slope

> v:= v+h*k: #improved Euler update for v

> t:= t+h: #update t
if frac(i/10)=0

> then print(t,v,sol(t),sol(t)-v);
fi;

> od:

>

So both the n=50 and n=100 improved Euler agree to 2 decimal places, over the entire interval.

One way to get better data about when the bolt reaches maximum height (v=0), and its contact velocity at time t=9.41 seconds would be to rerun improved Euler, with n=100, but only over the time intervals t=4..5 and t=0..10. One can use initial values from the table:

> t:=4.0;v:=5.530381164;n:=100;h:=1.0/n; #100 subdivisions,
#from t=4 to t=5.

> for i from 1 to n do

> k1:=f(t,v): #left-hand slope

> k2:=f(t+h,v+h*k1): #approximation to right-hand slope

> k:= (k1+k2)/2: #approximation to average slope

> v:= v+h*k: #improved Euler update for v

> t:= t+h: #update t
if t>4.5 and t<4.6 #I don't want all 100 entries

> then print(t,v,sol(t),sol(t)-v);
#also print the difference between exact solution
#and approximate solution.
fi;

> od:

This table makes it likely that the time of zero velocity it t=4.56 seconds, since the speed is smallest at that time.

Now for the velocity at time of impact:

>

> t:=9.0;v:=-39.88296265;n:=100;h:=1.0/n; #100 subdivisions,
#from t=9 to t=10.

> for i from 1 to n do

> k1:=f(t,v): #left-hand slope

> k2:=f(t+h,v+h*k1): #approximation to right-hand slope

> k:= (k1+k2)/2: #approximation to average slope

> v:= v+h*k: #improved Euler update for v

> t:= t+h: #update t
if t>9.35 and t<9.45 #I don't want all 100 entries

> then print(t,v,sol(t),sol(t)-v);
#also print the difference between exact solution
#and approximate solution.
fi;

> od:

So the velocity at time t=9.41 seconds is about -43.22 ft/sec. Of course, to find the time of impact without using the actual solution, we would have to numerically integrate v(t) to get y(t), and then find out when y(t) equaled zero.