VPF::[General] Нахождение производной - Форум программистов

FFFU

Дата 16.12.2010, 00:06 (ссылка)

(нет голосов)

Загрузка ...

Новичок

Профиль
Группа: Участник
Сообщений: 12
Регистрация: 23.10.2010

Репутация: нет
Всего: нет

Здравствуйте!
Нам в институте задали сделать программу, которая находит производные первого и второго порядка (приближенно). Я поискал в интернете, нашел такой хороший сайт.

Сделал все по аналогии, вот что получилось:

Код


      double precision nq, nq2, x, h
      double precision, parameter ::  pi=3.14159265358979
      integer i, n
      x=2
      h=1
      print*, "Сколько шагов?"
      read*, n
      print '(a10, a15, a30)', 'h', 'nq', 'nq2'

      do i=1, n
        nq=(f(x+2*h)-f(x+h))/h
        nq2=(f(x)-2*f(x+h)+f(x+2*h))/h**2
        print '(e10.4, f15.8, f30.8)', h, nq, nq2
        h=h/10
      end do

      contains
        function f(x)
        double precision f,x
          f=sin((pi*(x**2))/2)
        end function f
      end

Первую производную я посчитал на калькуляторе и он мне выдал 6,24544294. Программа, в принципе, близко к нему подходит, но точность все же хромает. А вторая производная вообще что-то непонятное выдает:

Код


 Сколько шагов?
13
         h             nq                           nq2
0.1000E+01    -0.99999930                   -1.99999913
0.1000E+00     3.68162834                  -23.22573704
0.1000E-01     6.30092851                    0.61995272
0.1000E-02     6.28760801                    2.89310373
0.1000E-03     6.28365383                    3.11677623
0.1000E-04     6.28323258                    3.13910620
0.1000E-05     6.28319019                    3.13945576
0.1000E-06     6.28318597                    3.46386414
0.1000E-07     6.28318553                    8.88177520
0.1000E-08     6.28318553                 -888.17844617
0.1000E-09     6.28318730                    0.00000000
0.1000E-10     6.28324059              8881784.19700125
0.1000E-11     6.28386232                    0.00000000

Помогите, пожалуйста. В чем моя ошибка?

Это сообщение отредактировал(а) FFFU - 16.12.2010, 00:07

kemiisto

Дата 16.12.2010, 13:45 (ссылка)

(голосов:1)

Загрузка ...

Дикий Кот. =^.^=

Награды: 1

Профиль
Группа: Участник Клуба
Сообщений: 3292
Регистрация: 29.7.2007

Репутация: нет
Всего: 160

Цитата(FFFU @ 15.12.2010, 22:06

)

нашел такой хороший сайт

УГ это, а не "хороший сайт".

Разбирать алгоритм ленно. А вот в реализации обращает на себя внимание вот эта строчка

Цитата(FFFU @ 15.12.2010, 22:06

)

double precision, parameter :: pi=3.14159265358979

Ошибка. Не будет тут двойной точности. Попробуйте вот такую программульку

Код


PROGRAM main

  IMPLICIT NONE

  DOUBLE PRECISION, PARAMETER :: pi = 3.14159265358979
  DOUBLE PRECISION, PARAMETER :: another_pi = 3.14159265358979D0

  PRINT *, pi
  PRINT *, another_pi

END PROGRAM main

и многое станет ясно.

Можно ещё так:

Код


  INTEGER, PARAMETER :: double = kind(0.0D0)
  REAL(kind=double), PARAMETER :: double_pi = 3.14159265358979D0

Добавлено через 1 минуту и 21 секунду
Или так:

Код


  INTEGER, PARAMETER :: double = kind(0.0D0)
  REAL(kind=double), PARAMETER :: double_pi = 3.14159265358979_double

Последний вариант, пожалуй, самый выгодный в долгосрочной перспективе.

--------------------

XMPP-канал.

Unit testing isn't enough. You need static typing.

Пей, кури, пиши на си! smile

FFFU

Дата 16.12.2010, 23:26 (ссылка)

(нет голосов)

Загрузка ...

Новичок

Профиль
Группа: Участник
Сообщений: 12
Регистрация: 23.10.2010

Репутация: нет
Всего: нет

Спасибо за совет, поменял.

Может быть и сайт подскажете, где можно посмотреть хороший пример?

0 Пользователей читают эту тему (0 Гостей и 0 Скрытых Пользователей)
0 Пользователей:
« Предыдущая тема \| Fortran \| Следующая тема »