Groups | Search | Server Info | Keyboard shortcuts | Login | Register [http] [https] [nntp] [nntps]
Groups > comp.lang.forth > #20489
| From | mhx@iae.nl (Marcel Hendrix) |
|---|---|
| Subject | Re: pde2 floating point benchmark code |
| Newsgroups | comp.lang.forth |
| Message-ID | <19669094008434@frunobulax.edu> (permalink) |
| Date | 2013-03-09 13:37 +0200 |
| References | <5a243597-c8bf-4e7c-8611-a5a88245429c@l16g2000yqe.googlegroups.com> |
| Organization | Wanadoo |
Krishna Myneni <krishna.myneni@ccreweb.org> writes Re: pde2 floating point benchmark code
> On Mar 8, 1:04 pm, m...@iae.nl (Marcel Hendrix) wrote:
>> Krishna Myneni <krishna.myn...@ccreweb.org> writes Re: pde2 floating point benchmark code
[..]
>>> However, speed is only one concern with this method of calculation.
>>> The grid size, set by dt and dx, must be small to obtain reasonable
>>> accuracy of the solution, even for short time intervals such as t = 1.
>>
>> Yes, that is really a problem. What is the error of this method?
>>
> Don't know yet. I hope to investigate this presently. It's a good
> method to illustrate the use of computing for solving partial
> differential equations, since the method is very easy to understand
> (just replace the partial derivatives with finite difference
> approximations).
Here's an experimental approach.
The outcome on my system is that the *difference between iterations*
goes down by a factor of four when dx is halved and dt adjusted so that
fx = 0.5. Each of these iteration takes eight times longer than the
previous one.
This does not say anything about the absolute / relative accuracy, but
allows one to stop when the solver makes too little progress. E.g.
when u_ext is 20 degrees C, a good stop would be when the difference
between iterations is less than 0.01 degree.
FORTH> 1e-5 GO
dx = 1.000000e-1 dt = 1.000000e-2 error 4.153411e0 0.001 seconds elapsed.
dx = 5.000000e-2 dt = 2.500000e-3 error 6.851562e-2 0.002 seconds elapsed.
dx = 2.500000e-2 dt = 6.250000e-4 error 1.693734e-2 0.004 seconds elapsed.
dx = 1.250000e-2 dt = 1.562500e-4 error 4.222571e-3 0.023 seconds elapsed.
dx = 6.250000e-3 dt = 3.906250e-5 error 1.054911e-3 0.168 seconds elapsed.
dx = 3.125000e-3 dt = 9.765625e-6 error 2.636820e-4 1.349 seconds elapsed.
dx = 1.562500e-3 dt = 2.441406e-6 error 6.591764e-5 10.810 seconds elapsed.
dx = 7.812500e-4 dt = 6.103516e-7 error 1.647923e-5 86.704 seconds elapsed.
dx = 3.906250e-4 dt = 1.525879e-7 error 4.119797e-6 696.044 seconds elapsed. ok
-marcel
-- ---------
ANEW -diffuse2
0.01e FVALUE dx \ 0.01 m
0.5e FVALUE D \ Diffusion coefficient
20e FVALUE u_ext
0.5e FVALUE fk \ choose dt so that fk < 0.5
dx FSQR f2/ D F/ FVALUE dt
1e fk F2* F- FVALUE 1-2fk
#10 VALUE nx
1e-5 FCONSTANT dxmin
10e dxmin F/ F>S 1+ =: nxmax
CREATE x nxmax DFLOATS ALLOT
CREATE u_x_0 nxmax DFLOATS ALLOT
CREATE u_j nxmax DFLOATS ALLOT
CREATE u_jp1 nxmax DFLOATS ALLOT
CREATE golden 0e DF, 0e DF, 0e DF,
\ Compute the solution at t=1
\ Compute the absolute difference with the previous iteration
: res? ( F: -- err )
CR ." dx = " dx E. ." dt = " dt E. ." error "
golden 0 DFLOAT[] DF@ u_j nx 7 #10 */ DFLOAT[] DF@ FDUP golden 0 DFLOAT[] DF! F- FSQR
golden 1 DFLOAT[] DF@ u_j nx 8 #10 */ DFLOAT[] DF@ FDUP golden 1 DFLOAT[] DF! F- FSQR F+
golden 2 DFLOAT[] DF@ u_j nx 9 #10 */ DFLOAT[] DF@ FDUP golden 2 DFLOAT[] DF! F- FSQR F+
3e F/ FSQRT FDUP E. SPACE ;
: evolve ( n -- ) ( F: -- res )
2/
0 ?DO
nx 1- 1 DO
u_j I 1+ DFLOAT[] DF@
u_j I 1- DFLOAT[] DF@ F+ fk F*
u_j I DFLOAT[] DF@ 1-2fk F* F+
u_jp1 I DFLOAT[] DF!
LOOP
u_jp1 1 DFLOAT[] DF@ u_jp1 0 DFLOAT[] DF! \ b.c.: u_x(0, t) = 0
u_ext u_jp1 nx 1- DFLOAT[] DF! \ b.c.: u(10, t) = u_ext
nx 1- 1 DO
u_jp1 I 1+ DFLOAT[] DF@
u_jp1 I 1- DFLOAT[] DF@ F+ fk F*
u_jp1 I DFLOAT[] DF@ 1-2fk F* F+
u_j I DFLOAT[] DF!
LOOP
u_j 1 DFLOAT[] DF@ u_j 0 DFLOAT[] DF! \ b.c.: u_x(0, t) = 0
u_ext u_j nx 1- DFLOAT[] DF! \ b.c.: u(10, t) = u_ext
LOOP res? ;
: INIT ( F: dx -- )
TO dx
dx FSQR F2/ D F/ TO dt
10e dx F/ F>S 1+ TO nx
nx nxmax > ABORT" init :: dx too small for preallocated array size"
0e nx 0 DO FDUP x I DFLOAT[] DF! dx F+ LOOP FDROP
nx 0 DO 1e u_x_0 I DFLOAT[] DF! LOOP
u_x_0 u_j nx DFLOATS MOVE ;
-- Stop when difference between two last iterations smaller than ERROR
: GO ( F: error -- )
FLOCAL error
0.1e FLOCAL xstart
golden 3 DFLOATS ERASE
BEGIN
xstart INIT
TIMER-RESET 1e ( second) dt F/ F>S EVOLVE .ELAPSED
error F> EKEY? 0= AND
WHILE
xstart F2/ TO xstart
REPEAT ;
0.01e init CR TIMER-RESET 1e dt F/ F>S evolve FDROP .ELAPSED
Back to comp.lang.forth | Previous | Next — Previous in thread | Next in thread | Find similar | Unroll thread
pde2 floating point benchmark code krishna.myneni@ccreweb.org - 2013-03-04 18:02 -0800
Re: pde2 floating point benchmark code Mark Wills <markrobertwills@yahoo.co.uk> - 2013-03-04 23:23 -0800
Re: pde2 floating point benchmark code krishna.myneni@ccreweb.org - 2013-03-05 05:22 -0800
Re: pde2 floating point benchmark code Doug Hoffman <glidedog@gmail.com> - 2013-03-05 12:11 -0500
Re: pde2 floating point benchmark code krishna.myneni@ccreweb.org - 2013-03-05 16:11 -0800
Re: pde2 floating point benchmark code Doug Hoffman <glidedog@gmail.com> - 2013-03-06 04:26 -0500
Re: pde2 floating point benchmark code krishna.myneni@ccreweb.org - 2013-03-06 05:22 -0800
Re: pde2 floating point benchmark code Doug Hoffman <glidedog@gmail.com> - 2013-03-06 09:44 -0500
Re: pde2 floating point benchmark code krishna.myneni@ccreweb.org - 2013-03-05 16:17 -0800
Re: pde2 floating point benchmark code mhx@iae.nl (Marcel Hendrix) - 2013-03-06 23:06 +0200
Re: pde2 floating point benchmark code krishna.myneni@ccreweb.org - 2013-03-07 04:57 -0800
Re: pde2 floating point benchmark code m.a.m.hendrix@tue.nl - 2013-03-07 05:49 -0800
Re: pde2 floating point benchmark code mhx@iae.nl (Marcel Hendrix) - 2013-03-07 21:12 +0200
Re: pde2 floating point benchmark code krishna.myneni@ccreweb.org - 2013-03-07 16:07 -0800
Re: pde2 floating point benchmark code Doug Hoffman <glidedog@gmail.com> - 2013-03-10 09:17 -0400
Re: pde2 floating point benchmark code mhx@iae.nl (Marcel Hendrix) - 2013-03-10 18:33 +0200
Re: pde2 floating point benchmark code Doug Hoffman <glidedog@gmail.com> - 2013-03-11 04:54 -0400
Re: pde2 floating point benchmark code humptydumpty <ouatubi@gmail.com> - 2013-03-07 22:32 -0800
Re: pde2 floating point benchmark code Krishna Myneni <krishna.myneni@ccreweb.org> - 2013-03-08 05:19 -0800
Re: pde2 floating point benchmark code mhx@iae.nl (Marcel Hendrix) - 2013-03-08 21:04 +0200
Re: pde2 floating point benchmark code Krishna Myneni <krishna.myneni@ccreweb.org> - 2013-03-08 16:31 -0800
Re: pde2 floating point benchmark code Zbiggy <zbigniew2011REMOVE@gmail.REMOVE.com> - 2013-03-09 00:56 +0000
Re: pde2 floating point benchmark code Krishna Myneni <krishna.myneni@ccreweb.org> - 2013-03-08 19:52 -0800
Re: pde2 floating point benchmark code Zbiggy <zbigniew2011REMOVE@gmail.REMOVE.com> - 2013-03-09 09:09 +0000
Re: pde2 floating point benchmark code mhx@iae.nl (Marcel Hendrix) - 2013-03-09 13:37 +0200
Re: pde2 floating point benchmark code Krishna Myneni <krishna.myneni@ccreweb.org> - 2013-03-09 06:28 -0800
Re: pde2 floating point benchmark code mhx@iae.nl (Marcel Hendrix) - 2013-03-09 17:58 +0200
Re: pde2 floating point benchmark code humptydumpty <ouatubi@gmail.com> - 2013-03-08 12:23 -0800
Re: pde2 floating point benchmark code Krishna Myneni <krishna.myneni@ccreweb.org> - 2013-03-08 16:43 -0800
csiph-web