!********************************************************* !* Calculate a limited development of a real function * !* f(x)/g(x) at x=0 up to order 25, knowing the limited * !* developments of f(x) and g(x). * !* ----------------------------------------------------- * !* SAMPLE RUN: * !* * !* Limited development of a real function f(x)/g(x): * !* * !* Function to develop: exp(x)/1+x * !* * !* Limited development of f(x)=exp(x) is: * !* 1 + x + x^2/2! + x^3/3! + ... + x^n/n! + ... * !* * !* Limited development of g(x)=1+x is: * !* 1 + x * !* * !* Order of development (max=25): 6 * !* * !* Input coefficients of limited dev. of f(x): * !* a0 = 1 * !* a1 = 1 * !* a2 = 0.5 * !* a3 = 0.1666667 * !* a4 = 0.0416667 * !* a5 = 0.0083333 * !* a6 = 0.0013889 * !* * !* Input coefficients of limited dev. of g(x): * !* b0 = 1 * !* b1 = 1 * !* b2 = 0 * !* b3 = 0 * !* b4 = 0 * !* b5 = 0 * !* b6 = 0 * !* * !* The coefficients of limited dev. of f(x)/g(x) are: * !* c0 = 1.0000000 * !* c1 = 0.0000000 * !* c2 = 0.5000000 * !* c3 = -0.3333333 * !* c4 = 0.3750000 * !* c5 = -0.3666667 * !* c6 = 0.3680556 * !* * !* Function to develop: 1/1+sh(x) * !* * !* Limited development of g(x)=1+sh(x) is: * !* 1 + x + x^3/3! + x^5/5! + ... + x^2n+1/(2n+1)! + ... * !* * !* Order of development (max=25): 6 * !* * !* Input coefficients of limited dev. of f(x): * !* a0 = 1 * !* a1 = 0 * !* a2 = 0 * !* a3 = 0 * !* a4 = 0 * !* a5 = 0 * !* a6 = 0 * !* * !* Input coefficients of limited dev. of g(x): * !* b0 = 1 * !* b1 = 1 * !* b2 = 0 * !* b3 = 0.1666667 * !* b4 = 0 * !* b5 = 0.0083333 * !* b6 = 0 * !* * !* The coefficients of limited dev. of f(x)/g(x) are: * !* c0 = 1.0000000 * !* c1 = -1.0000000 * !* c2 = 1.0000000 * !* c3 = -1.1666667 * !* c4 = 1.3333334 * !* c5 = -1.5083334 * !* c6 = 1.7111112 * !* * !* ----------------------------------------------------- * !* Ref.: "Mathematiques en Turbo Pascal by Alain * !* Reverchon and Marc Ducamp, Editions Eyrolles, * !* Paris, 1991" [BIBLI 03]. * !* * !* F90 version by J-P Moreau. * !* (www.jpmoreau.fr) * !********************************************************* Program LTDDEV2 parameter(NSIZE=25) integer i,m real*8 T1(0:NSIZE),T2(0:NSIZE),R(0:NSIZE) print *,' ' print *,' Limited development of a real function f(x)/g(x):' print *,' ' print *,' Function to develop: exp(x)/1+x' print *,' ' print *,' Limited development of f(x)=exp(x) is:' print *,' 1 + x + x^2/2! + x^3/3! + ... + x^n/n! + ...' print *,' ' print *,' Limited development of g(x)=1+x is:' print *,' 1 + x' print *,' ' write(*,10,advance='no'); read *, m print *,' ' print *,' Input coefficients of limited dev. of f(x):' do i=0, m write(*,20,advance='no') i; read *, T1(i) end do print *,' ' print *,' Input coefficients of limited dev. of g(x):' do i=0, m write(*,21,advance='no') i; read *, T2(i) end do print *,' ' call ar_dldiv(m,T1,T2,R) print *,' The coefficients of limited dev. of f(x)/g(x) are:' do i=0, m write(*,30) i, R(i) end do print *,' ' 10 format(' Order of development (max=25): ') 20 format(' a',I1,' = ') 21 format(' b',I1,' = ') 30 format(' c',I1,' = ',F10.7) END Subroutine ar_dldiv(n, t1, t2, res) parameter(NSIZE=25) real*8 x,t1(0:NSIZE),t2(0:NSIZE),res(0:NSIZE) if (n > NSIZE) return if (dabs(t2(0)) < 1.d-16) return res(0)=t1(0)/t2(0) do i=1, n x = t1(i) do j=1, i x = x-t2(j)*res(i-j) end do res(i) = x end do return End !end of file ltddev2.f90