Groups | Search | Server Info | Keyboard shortcuts | Login | Register [http] [https] [nntp] [nntps]
Groups > comp.lang.forth > #11962
| From | mhx@iae.nl (Marcel Hendrix) |
|---|---|
| Subject | Re: code for reading gif/png files |
| Newsgroups | comp.lang.forth |
| Message-ID | <80141496988435@frunobulax.edu> (permalink) |
| Date | 2012-05-07 21:21 +0200 |
| References | <b60188e4-ff82-44b7-ba73-d94af2aaf61d@p6g2000yqi.googlegroups.com> |
| Organization | Wanadoo |
Krishna Myneni <krishna.myneni@ccreweb.org> writes Re: FSL contribution available for review
> On Mar 6, 5:07 pm, m...@iae.nl (Marcel Hendrix) wrote:
>> m...@iae.nl (Marcel Hendrix) wrote Re: FSL contribution available for review
>> > Krishna Myneni <krishna.myn...@ccreweb.org> writes Re: FSL contribution available for review
[..]
> One of the nice properties of the Numerov difference equation is that
> the error in the solution goes as O(h^4), where h is the step size
> (see, for example, http://research.physics.illinois.edu/ElectronicStructure/498CQM/lnotes/lec2A.html).
> Marcel's error calculation with the MPFR library directly shows the
> h^4 error dependence.
[..]
With a small modification Numerov becomes O(h^8).
\ The maximum residual from R = 1.000000e-3 to R = 1.000000e1 is 5.101759e-8 ( n=250 ) ok
\ The maximum residual from R = 1.000000e-3 to R = 1.000000e1 is 1.981971e-10 ( n=500 ) ok
\ The maximum residual from R = 1.000000e-3 to R = 1.000000e1 is 7.837742e-13 ( n=1,000 ) ok
\ The maximum residual from R = 1.000000e-3 to R = 1.000000e1 is 1.463233e-14 ( n=2,000 ) ok
-marcel
-- -------------------------------------------------------------
NEEDS -fsl_util
ANEW -numerov8
FALSE
[IF]
Integration of the 2nd order differential equation
y''(x) = f(x,y)
with Q(x) a given function, using the Numerov algorithm.
The Numerov algorithm may be expressed as shown in (1),
and then has an error varying like the eight power of h.
This implementation uses a function which evaluates f(x,y),
the initial and final argument values, and the number of steps
for the integration. The stepsize h is computed.
Usage:
Define a word myF ( f: x y -- F{x,y} )
' myF IS Ffunc
Prepare arrays, using the array size and the integration limits
n Xstart Xend setup
Establish initial conditions and Y values
Y_0 Y_1 setY
Perform integration
NumInt
Results are in the array Y
[1] 'Explicit eighth order two-step methods with nine
stages for integrating oscillatory problems', Ch. Tsitouras.
[THEN]
#10 =: st \ stages
st FLOAT ARRAY c{
c{ st }FREAD
-1e 0e -1.618033988749895e -0.08935969452190693e -0.7180027509073757e 0.7180027509073757e -0.25e 0.25e -1e 1e
st FLOAT ARRAY b{
b{ st }FREAD
0.02267478608411768e 0e 0e 0e 0.1091598371161353e 0.1091598371161353e 0.3880338950775969e 0.3880338950775969e -0.01986851827784987e 0.002806267806267806e
st st FLOAT MATRIX a{{
a{{ st st }}FREAD
0e 0e 0e 0e 0e 0e 0e 0e 0e 0e
0e 0e 0e 0e 0e 0e 0e 0e 0e 0e
0.4363389981249825e 0.06366100187501753e 0e 0e 0e 0e 0e 0e 0e 0e
-0.02663944838475621e -0.02138085097354293e 0.007333029599869930e 0e 0e 0e 0e 0e 0e 0e
-0.05259994463359025e 0.1179873479656171e 0.006223764486158627e -0.1728485681165938e 0e 0e 0e 0e 0e 0e
-0.1594931414841811e 1.756644381705087e 0.002177668974400012e -1.462560200318788e 0.4799966417324492e 0e 0e 0e 0e 0e
-0.01315251843525407e 0.08148753879227717e 0.002255441346558031e -0.1407999204529257e -0.02359301393743279e 0.00005247268677732879e 0e 0e 0e 0e
0.1182251406950030e -0.2071467658425108e -0.009902612273876664e 0.2377506314405291e -0.1720715921748083e 0.008456715906120000e 0.1809384822495436e 0e 0e 0e
0.6545342597532786e 4.968502507588174e -0.05384950599580273e -4.016696408666935e -1.055358930155700e 0.2067362330539400e 1.043495190976432e -1.747363346553386e 0e 0e
-0.2731258141928670e -19.26209659195308e 0.2868033393908071e 21.50877058850632e -1.286133152186278e 0.7520725477949123e -1.229894203564763e 0.6765130737370460e -0.1729097875320912e 0e
\ Initialization
0 VALUE n ( 1 + number of steps / integer greater than 1 )
FALSE VALUE used? ( have arrays already been allocated? )
0e FVALUE h ( stepsize )
FLOAT DARRAY Y{
FLOAT DARRAY X{
st FLOAT ARRAY F{
DEFER Ffunc
:NONAME ( f: x y -- f{x} ) F2DROP 0e ; IS Ffunc
\ Allocate arrays and compute h and X values
: setarrays ( n -- )
TO n
used? IF X{ }free Y{ }free ENDIF
X{ n }malloc Y{ n }malloc
TRUE TO used? ;
: stepsize ( f: Xstart Xend -- Xstart h ) FOVER F- n 1- S>F F/ FDUP TO h ;
: fillX ( f: Xstart h -- )
X{ n 0 DO DUP F2DUP ( f: Xstart h Xstart h )
I S>F F* F+ I } F! ( f: Xstart h )
LOOP DROP F2DROP ;
: F*b ( F: -- r ) 0e st 0 ?DO F{ I } F@ b{ I } F@ F* F+ LOOP ;
: F*a(o) ( row -- ) ( F: -- r ) LOCALS| row | 0e st 0 ?DO F{ I } F@ a{{ row I }} F@ F* F+ LOOP ;
: setup ( n -- ) ( f: Xstart Xend -- ) setarrays stepsize fillX ;
: setY ( f: Y_0 Y_1 -- ) Y{ 1 } F! Y{ 0 } F! ;
: NumInt ( -- )
0e FLOCAL f0
X{ 0 } F@ Y{ 0 } F@ Ffunc FLOCAL f1
h FSQR FLOCAL h^2
n 1- 1
?DO \ first stage
f1 TO f0 f0 F{ 0 } F!
X{ I } F@ Y{ I } F@ Ffunc TO f1 f1 F{ 1 } F!
st 2 ?DO \ another 8 stages
c{ I } F@ h F* X{ J } F@ F+
c{ I } F@ F1+ Y{ J } F@ F*
c{ I } F@ Y{ J 1- } F@ F* F-
h^2 I F*a(o) F* F+
( x y ) Ffunc F{ I } F!
LOOP
Y{ I } F@ F2* Y{ I 1- } F@ F- h^2 F*b F* F+ Y{ I 1+ } F!
LOOP ;
Back to comp.lang.forth | Previous | Next — Previous in thread | Next in thread | Find similar | Unroll thread
code for reading gif/png files Krishna Myneni <krishna.myneni@ccreweb.org> - 2012-05-07 03:56 -0700
Re: code for reading gif/png files mhx@iae.nl (Marcel Hendrix) - 2012-05-07 21:21 +0200
Re: code for reading gif/png files mhx@iae.nl (Marcel Hendrix) - 2012-05-07 21:27 +0200
Re: code for reading gif/png files Hans Bezemer <the.beez.speaks@gmail.com> - 2012-05-07 23:36 +0200
Re: code for reading gif/png files Krishna Myneni <krishna.myneni@ccreweb.org> - 2012-05-07 15:12 -0700
Re: code for reading gif/png files Krishna Myneni <krishna.myneni@ccreweb.org> - 2012-05-07 16:07 -0700
Re: code for reading gif/png files Tarkin <tarkin000@gmail.com> - 2012-05-07 17:31 -0700
Re: code for reading gif/png files Krishna Myneni <krishna.myneni@ccreweb.org> - 2012-05-07 18:28 -0700
Re: code for reading gif/png files Rugxulo <rugxulo@gmail.com> - 2012-05-07 21:15 -0700
Re: code for reading gif/png files Krishna Myneni <krishna.myneni@ccreweb.org> - 2012-05-09 05:25 -0700
Re: code for reading gif/png files Rugxulo <rugxulo@gmail.com> - 2012-05-09 15:17 -0700
Re: code for reading gif/png files Krishna Myneni <krishna.myneni@ccreweb.org> - 2012-05-11 03:33 -0700
csiph-web