>> From: Maple Group
-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-
Date: Thu, 26 Oct 2000 14:40:24 -0700 (PDT)
From: Robert Israel
To:
Subject: Runge Kutta order 4
On Wed, 25 Oct 2000, Chuck Baker wrote:
> Is it possible to solve the following using Maple 6?
Yes, it is. But you shouldn't be asking this group to do your
homework for you.
Robert Israel
Department of Mathematics http://www.math.ubc.ca/~israel
University of British Columbia
Vancouver, BC, Canada V6T 1Z2
-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-
Date: Thu, 26 Oct 2000 15:45:56 -0700
To:
From: David Harrington
Subject: Runge Kutta order 4
This will do it all in Maple, but there is no difference between the real
and numerical solutions at the 4 digit accuracy with the rk4 method
(substitute foreuler for rk4 to see a difference, or use a bigger stepsize).
The formatting could be made nicer by using the printf command.
restart;de:={diff(y(x),x)=1-x+4*y(x),y(0)=1};
d
de := {-- y(x) = 1 - x + 4 y(x), y(0) = 1}
dx
> yanal:=rhs(dsolve(de,y(x)));
19
yanal := - 3/16 + 1/4 x + -- exp(4 x)
16
> ynum:=dsolve(de,y(x),numeric,method=classical[rk4],'stepsize'=0.01);
ynum := proc(x_classical) ... end proc
> ynum(1);
[x = 1, y(x) = 64.8977978130418052]
> for xx from 0 to 1 by 1/10 do
> evalf(xx,1),evalf(subs(x=xx,yanal),5),evalf(rhs(ynum(xx)[2]),5);
> od;
0., 1.0000, 1.
.1, 1.6090, 1.6090
.2, 2.5053, 2.5053
.3, 3.8301, 3.8301
.4, 5.7942, 5.7942
.5, 8.7121, 8.7120
.6, 13.052, 13.053
.7, 19.516, 19.516
.8, 29.146, 29.145
.9, 43.498, 43.498
1., 64.898, 64.898
-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-
Date: Thu, 26 Oct 2000 23:07:52 -0400 (EDT)
From: Carl DeVore
To:
Subject: Runge Kutta order 4
Yes, the whole thing is trivial in Maple and can be written in 10-15
minutes.
> restart;
> with(linalg):
> Digits:= trunc(evalhf(Digits)):
> IVP:= {diff(y(x), x) = 1-x+4*y(x), y(0)=1}:
> X:= array(1..101, [seq(.01*k, k= 0..100)]):
> Y[analytical]:= map(unapply(rhs(dsolve(IVP, {y(x)})), x), X):
> Y[RK4]:= col(dsolve
> (IVP
> ,{y(x)}
> ,type= numeric
> ,method= classical[rk4]
> ,value= X
> ,stepsize= .01
> )[2,1]
> ,2):
> stackmatrix
> (convert([`x-value`, `y-analytical`, `y-RK4`, `Abs. Error`], list)
> ,augment
> (X
> ,evalf(eval(Y[analytical]), 6)
> ,evalf(eval(Y[RK4]), 6)
> ,evalf(map(abs, evalm(Y[analytical]-Y[RK4])), 10)
> )
> );
... output deleted
Carl Devore
University of Delaware
-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-
Date: Fri, 27 Oct 2000 07:56:13 -0400 (EDT)
From: John Little
Subject: Runge Kutta order 4
To:
Dear Chuck,
Don't mean to brag for Maple, but this kind of thing is
almost too easy -- a perfect illustration of the way Maple can
combine symbolic and numerical computation. You could also
plot the solution of your ODE in various ways.
Here's a Maple V procedure that generates the table of
exact solution values, Runge-Kutta approximations, and absolute
errors you asked for. I did this more or less ``off the top of my
head'', so some Maple purists might find better ways to do some of
the steps. I also didn't try to make the output pretty
and I didn't worry about trying to print out only the first
four decimal places in the numbers. About the numerical aspects:
Maple's default floating point number system uses 10 decimal digits
and the output was generated that way. There is also a
system variable called Digits you can set if you want to use
more or fewer decimal digits in floating point numbers.
RungeKuttaEx := proc()
#
# Generate approximate solution of y' = 1 - x + 4*y; y(0) = 1
# by the Classical RungeKutta Method and compare with analytic
# solution
#
#
local i, h, k1, k2, k3, k4 ,t, w, lastt, lastw, ansol, anval, f;
# first find the analytic solution of the IVP using dsolve
ansol := subs(dsolve({diff(y(x),x)=1+x-4*y(x),y(0)=1},y(x)),y(x));
# now set up for the Runge-Kutta method
f := (x,y) -> 1 + x - 4*y;
h := .01;
lastt := 0;
lastw := 1;
for i to 100 do
t := lastt + h;
k1 := f(lastt,lastw)*h;
k2 := f(lastt + h/2,lastw + k1/2)*h;
k3 := f(lastt + h/2,lastw + k2/2)*h;
k4 := f(t,lastw + k3)*h;
w := lastw + (k1 + 2*k2 + 2*k3 + k4)/6;
if evalb(i mod 10 = 0) then
anval:=evalf(subs(x=t,ansol));
lprint(`t=`,t,`exact value=`,anval,
`R-K approximation=`,w,`error=`,abs(w-anval));
fi;
lastw := w;
lastt := t:
od:
RETURN()
end:
The output:
> RungeKuttaEx();
t= .10 exact value= .7571350374 R-K approximation= .7571350422
error= .48e-8
t= .20 exact value= .6025797833 R-K approximation= .6025797899
error= .66e-8
t= .30 exact value= .5072202972 R-K approximation= .5072203037
error= .65e-8
t= .40 exact value= .4515409209 R-K approximation= .4515409265
error= .56e-8
t= .50 exact value= .4224599176 R-K approximation= .4224599226
error= .50e-8
t= .60 exact value= .4112083371 R-K approximation= .4112083411
error= .40e-8
t= .70 exact value= .4119081759 R-K approximation= .4119081790
error= .31e-8
t= .80 exact value= .4206192907 R-K approximation= .4206192930
error= .23e-8
t= .90 exact value= .4347005245 R-K approximation= .4347005261
error= .16e-8
t= 1.00 exact value= .4523814566 R-K approximation= .4523814579
error= .13e-8
John B. Little
Dept. of Mathematics and Computer Science
College of the Holy Cross
Worcester, MA 01610 USA
(508) 793-2274
email:
homepage: http://math.holycross.edu/~little/homepage.html
|
blogs, forums etc