+ Ответить в теме
Показано с 1 по 5 из 5

Тема: Интерполяционные многочлены Лагранжа - Fortran

  1. #1
    dummy Turb_man is on a distinguished road
    Регистрация
    30.06.2017
    Возраст
    22
    Сообщений
    9
    Вес репутации
    0

    По умолчанию Интерполяционные многочлены Лагранжа - Fortran

    Здравствуйте уважаемые пользователи форума.
    У меня возникла одна небольшая проблема.
    Передо мной стояла следующая задача: найти зависимость скорости частицы(или любого другого объекта) в произвольной точке.
    Значение в узлах сетки я считываю из предварительно данных мне файлов.
    То есть у меня есть функция в трёх переменных в трёхмерном пространстве, которую нужно интерполировать с помощью полиномов Лагранжа
    по заданным значениям в узлах сетки.
    Случай интерполяции функции одной переменной не вызвал у меня затруднений. А вот здесь я немного сел в лужу.

    http://stu.sernam.ru/book_dig_m.php?id=15

    Это ссылка на формулу многочленов Лагранжа в двухмерном случае.

    http://pmpu.ru/vf4/interpolation

    А здесь случай произвольной размерности.

    В трёхмерном случае узлы интерполяции, не должны лежать в одной плоскости. И здесь я поэтому запутался немного с циклами. Если кто-то найдёт ошибку или уже сталкивался с подобной проблей, то сообщите пожалуйста.

    Уже месяц периодически мучаюсь с этой программой, которую писал сам, так что сразу извиняюсь за корявость кода.

    Немного поясню суть метода Лагрнажа.

    Пусть у нас есть область. В данном случае это прямоугольный параллелепипед. Мы разбиваем его на одинаковые маленькие прямоугольные паралеллепипеды. То есть получается такая трёхмерная сетка. В каждом узле сетки заданы значение функции. Теперь с помощью интерполяционных полиномов Лагранжа следует найти значение функции в произвольной точке.

    Привожу код моей программы. Программа написана на языке Fortran.

  2. По умолчанию

     
    Хотите избавиться от рекламы? Зарегистрируйтесь
  3. #2
    dummy Turb_man is on a distinguished road
    Регистрация
    30.06.2017
    Возраст
    22
    Сообщений
    9
    Вес репутации
    0

    По умолчанию Re: Интерполяционные многочлены Лагранжа - Fortran

    Код fortran:
    1. program time
    2.  
    3.       call v_time
    4.  
    5.       end
    6. c------------------------------------------------------------------------
    7. c Определение  номера элемента, в который попадёт точка  
    8.       subroutine place(xx,yy,zz,xstep,ystep,zstep)
    9.       real*8 zz(1), xx(1), yy(1), step, h, w, s
    10.       integer m, n, k, l, xstep(1), ystep(1), zstep(1)
    11. c Задание точек, в которых нужно найти скорость жидкости
    12.       step = 0.1
    13. c Размеры исходной области
    14.       h = 4 * acos(-1.0)/3
    15.       w = 1
    16.       s = 4 * acos(-1.0)
    17.       do m = 1, 1
    18.         xx(1) = 0.0
    19.         yy(1) = 0.0
    20.         zz(1) = 0.0
    21. c        xx(m + 1) = xx(m) + step
    22. c        yy(m + 1) = yy(m) + step
    23. c        zz(m + 1) = zz(m) + step
    24.       enddo
    25. c Определение номера элемента по оси z, в который попадёт точка
    26.       do m = 1, 1
    27.         do n = 0, 24
    28.            if(zz(m).ge.(n*h/25) .and. (zz(m).le.((n + 1)*h/25)))
    29.      & then
    30.              zstep(m) = n + 1
    31.            endif
    32.         enddo
    33.       enddo
    34. c Определение номера элемента по оси x, в который попадёт точка
    35.       do m = 1, 1
    36.         do l = 0, 29
    37.            if(xx(m).ge.(l*s/30) .and. (xx(m).le.((l + 1)*s/30)))
    38.      & then
    39.              xstep(m) = l + 1
    40.            endif
    41.         enddo
    42.       enddo
    43. c Определение номера элемента по оси y, в который попадёт точка
    44.        do m = 1, 1
    45.         do k = 0, 9
    46.            if(yy(m).ge.(k*w/10) .and. (yy(m).le.((k + 1)*w/10)))
    47.      & then
    48.              ystep(m) = k + 1
    49.            endif
    50.        enddo
    51.        enddo
    52.       end
    53. c--------------------------------------------------------------------------
    54. c Вычисление базисных интерполяционных полиномов Лагранжа
    55.       subroutine lagrange(xx, yy, zz, x, y, z, l, i, j, k)
    56.       integer i, j, k, s, q, r
    57.       real*8  zz, xx, yy, x(8,8,8), y(8,8,8), z(8,8,8), l(8,8,8)
    58.       l = 1
    59.       do r = 1, 8
    60.         do s = 1, 8
    61.           do q = 1, 8
    62.             if(i.ne.r .and. j.ne.s .and. k.ne.q) then
    63.               l(i, j, k) = l(i, j, k) *
    64.      & * (xx - x(r,s,q))/(x(i,j,k) - x(r,s,q))
    65.      & * (yy - y(r,s,q))/(y(i,j,k) - y(r,s,q))
    66.      & * (zz - z(r,s,q))/(z(i,j,k) - z(r,s,q))
    67.             else
    68.                l(i, j, k) = l(i, j, k) * 1
    69.             endif
    70.  enddo
    71.         enddo
    72.       enddo
    73. c      write(*,*) l(i,j,k)
    74.       end
    75. c--------------------------------------------------------------------------
    76. c Сортировка точек в элементе с выделенной точкой
    77. c для применения базисных Полиномов
    78.       subroutine polpoint(pox, poy, poz, num_step)
    79.       integer xstep(1), ystep(1), zstep(1)
    80.       integer l, n, num_step(512)
    81.       integer r, q, s
    82.       integer size_file, n1, el, gel, ex, ey, ez, i, j, k
    83.       real*8 x, y, z, zz(1), xx(1), yy(1), t(512)
    84.       real*8 u(512), v(512), pox(1,8,8,8)
    85.       real*8 poy(1,8,8,8), poz(1,8,8,8)
    86.       character(10) second_name
    87. c Файл с расположением точек
    88.       second_name = 'Geom_total'
    89.  
    90.       call place(xx,yy,zz,xstep,ystep,zstep)
    91.  
    92.       open(unit = 1, file = second_name, form = 'unformatted',
    93.      & status='old')
    94.  
    95.       read(1) size_file
    96. c Я считываю данные из файла и записываю их так, чтобы применить метод Лагранжа.
    97.       do l = 1, size_file
    98.         read(1) n1, el, gel, ex, ey, ez, i, j, k, x,
    99.      &  y,  z
    100.         do n = 1, 1
    101.           if(ex.eq.xstep(n) .and. (ey-10).eq.ystep(n)) then
    102.             if(ez.eq.zstep(n)) then
    103.               do r = 1, 8
    104.                 if(i.eq.r) then
    105.                   do s = 1, 8
    106.                     if(j.eq.s) then
    107.                       do q = 1, 8
    108.                         if(k.eq.q) then
    109.                           t(q+8*(s-1)+64*(r-1)+512*(n-1)) = x
    110.                           u(q+8*(s-1)+64*(r-1)+512*(n-1)) = y
    111.                           v(q+8*(s-1)+64*(r-1)+512*(n-1)) = z
    112.                           num_step(q+8*(s-1)+64*(r-1)+512*(n-1)) = l
    113.                         endif
    114.                       enddo
    115.                     endif
    116.                   enddo
    117.                 endif
    118.               enddo
    119.             endif
    120.           endif
    121.         enddo
    122.       enddo
    123.       close(1)
    124. c Отсортированные точки в спекиральном элементе
    125.       do n = 1, 1
    126.         do r = 1, 8
    127.           do s = 1, 8
    128.             do q = 1, 8
    129.               pox(n,r,s,q) = t(q+8*(s-1)+64*(r-1)+512*(n-1))
    130.               poy(n,r,s,q) = u(q+8*(s-1)+64*(r-1)+512*(n-1))
    131.               poz(n,r,s,q) = v(q+8*(s-1)+64*(r-1)+512*(n-1))
    132. c              write(3,*) pox(m,r,s,q), poy(m,r,s,q), poz(m,r,s,q)
    133.             enddo
    134.           enddo
    135.         enddo
    136.       enddo
    137.       end
    Последний раз редактировалось AiK; 07.07.2017 в 13:38. Причина: [code=fortran][/code]

  4. #3
    dummy Turb_man is on a distinguished road
    Регистрация
    30.06.2017
    Возраст
    22
    Сообщений
    9
    Вес репутации
    0

    По умолчанию Re: Интерполяционные многочлены Лагранжа - Fortran

    c Сорртировка скоростей в узлах сетки, для применения базисных полиномов
    Код fortran:
    1.       subroutine velocity(num_step, third_name, v11, v22, v33)
    2.       character(16) third_name
    3.       real*8 vx, vy, vz, p
    4.       real*8 v1(512), v2(512), v3(512)
    5.       real*8 v11(1,8,8,8), v22(1,8,8,8), v33(1,8,8,8)
    6.       integer size_file, num_step(512), m, l, r, s, q
    7.       open(unit = 2, file = third_name, form='unformatted',
    8.      &  status='old')
    9.       read(2) size_file
    10. c Я считываю данные из файла и записываю их так, чтобы применить метод Лагранжа.
    11.       do l = 1, size_file
    12.           read(2) vx, vy, vz, p
    13.           do m = 1, 512
    14.             if(l.eq.num_step(m)) then
    15.                v1(m) = vx
    16.                v2(m) = vy
    17.                v3(m) = vz
    18.             endif
    19.           enddo
    20.       enddo
    21.  close(2)
    22. c Отсортированный массив трёх компонет скорости
    23.       do m = 1, 1
    24.         do r = 1, 8
    25.           do s = 1, 8
    26.             do q = 1, 8
    27.               v11(m,r,s,q) = v1(q+8*(s-1)+64*(r-1)+512*(m-1))
    28.               v22(m,r,s,q) = v2(q+8*(s-1)+64*(r-1)+512*(m-1))
    29.               v33(m,r,s,q) = v3(q+8*(s-1)+64*(r-1)+512*(m-1))
    30. c              write(3,*) v11(m,r,s,q), v22(m,r,s,q), v33(m,r,s,q)
    31.             enddo
    32.           enddo
    33.         enddo
    34.       enddo
    35.       end
    36. c--------------------------------------------------------------------------
    37. c--------------------------------------------------------------------------
    38.  Вычисление среднего значения скорости для каждой точки
    39.       subroutine av_velocity(xx, yy, zz, pox, poy, poz, v11, v22, v33)
    40.       real*8 v11(1,8,8,8), v22(1,8,8,8), v33(1,8,8,8), lag(1,8,8,8)
    41.       integer i, j, k, n, xstep(1), ystep(1), zstep(1)
    42.       real*8 vx(1), vy(1), vz(1), xx(1), yy(1), zz(1)
    43.       real*8 pox(1,8,8,8), poy(1,8,8,8), poz(1,8,8,8)
    44.       call place(xx, yy, zz, xstep, ystep, zstep)
    45.       vx = 0
    46.       vy = 0
    47.       vz = 0
    48.       lag = 1
    49.       do n = 1, 1
    50.         do i = 1, 8
    51.           do j = 1, 8
    52.             do k = 1, 8
    53. c              write(*,*) v11(n,i,j,k), v22(n,i,j,k), v33(n,i,j,k)
    54. c Применние интерполяционных полиномов Лагранжа
    55.               call lagrange(xx(n),yy(n),zz(n),pox(n,:,:,:),poy(n,:,:,:),
    56.      & poz(n,:,:,:), lag(n,:,:,:), i, j, k)
    57.               vx(n) = vx(n) + v11(n,i,j,k) * lag(n,i,j,k)
    58.               vy(n) = vy(n) + v22(n,i,j,k) * lag(n,i,j,k)
    59.               vz(n) = vz(n) + v33(n,i,j,k) * lag(n,i,j,k)
    60.             enddo
    61.           enddo
    62.         enddo
    63.       enddo
    64.       call place(xx, yy, zz, xstep, ystep, zstep)
    65.       do n = 1, 1
    66. c          write(3,*) xx(n), yy(n), zz(n), vx(n), vy(n), vz(n)
    67.       enddo
    68.       end
    69. c----------------------------------------------------------------------------
    70. c Нахождение зависимости эволюции средних скоростей
    71. c частиц в зависимости отвремени
    72.       subroutine v_time
    73.  
    74.       integer i, num_step(512)
    75.       character(6) first_name
    76.       character(16) third_name
    77.       real*8 v11(1,8,8,8), v22(1,8,8,8), v33(1,8,8,8), poz(1,8,8,8)
    78.       real*8 xx(1), yy(1), zz(1),  pox(1,8,8,8), poy(1,8,8,8)
    79. c      call polpoint(lag, num_step)
    80.  
    81.       third_name = 'Velo_total_00001'
    82.       first_name = '5_file'
    83.       open(unit = 3, file = first_name)
    84.       call polpoint(pox, poy, poz, num_step)
    85.       do i = 1, 1
    86.          write(third_name(12:16),'(I5.5)') i
    87.          call velocity(num_step, third_name, v11, v22, v33)
    88.          call av_velocity(xx, yy, zz, pox, poy, poz, v11, v22, v33)
    89.          write(3,*) " "
    90.       enddo
    91.       close(3)
    92.   end
    Последний раз редактировалось AiK; 07.07.2017 в 13:39. Причина: [code=fortran][/code]

  5. #4
    Administrator Админ
    system architect
    AiK is on a distinguished road Аватар для AiK
    Регистрация
    13.02.2004
    Адрес
    СПб
    Сообщений
    2,298
    Вес репутации
    80

    По умолчанию Re: Интерполяционные многочлены Лагранжа - Fortran

    найти зависимость скорости частицы(или любого другого объекта) в произвольной точке.
    Зависимость от чего?

    Немного поясню суть метода Лагрнажа.
    Немного описана суть интерполяции, про метод Лагранжа не сказано ни слова. Кроме того, скорость - это первая производная от пути. То есть не исключено применение интерполяционных полиномов Эрмита.

    Привожу код моей программы.
    Кода нет.
    Даже самый дурацкий замысел можно воплотить мастерски

  6. #5
    dummy Turb_man is on a distinguished road
    Регистрация
    30.06.2017
    Возраст
    22
    Сообщений
    9
    Вес репутации
    0

    По умолчанию Re: Интерполяционные многочлены Лагранжа - Fortran

    Спасибо за уточнения.

    1) Нужно найти скорость частицы в зависимости от времени.

    2) Тут я неправильно сформулировал, это суть интерполяции. Про метод Лагранжа можно почитать по ссылкам выше.

    3) У вас такие правила по ограничению количества символов, поэтому я и подгружал по частям.

+ Ответить в теме

Похожие темы

  1. а как же Fortran?
    :( Доброго времени суток.. Почему у Вас нет раздела для тех кто программирует на Фортране? У меня возникла проблема с графикой.. и я даже не знаю...
    от Spinyz в разделе Жалобная книга
  2. Связанные списки и многочлены; Delphi7
    Помогите составить алгоритм!! Никак не могу сообразить, как в линейном списке задать инфу о 2-х значениях!!:( Пусть многочлен задан в виде линейного...
    от _Kommandor_ в разделе задачи на Паскале и Delphi
  3. вопрос по Fortran
    Книга Бартеньева «Современный Фортран» является учебным пособием по языку программирования Microsoft Fortran PowerStation 4.0. И 12 глава его книги...
    от Spinyz в разделе Вопрошайка
  4. Fortran. Задача из курсовой
    Надеюсь не наделаю много шума) сам я в маи учусь, задача тут угребистая попалась, никак решить не могу. а 24ого кровь из носу нести надо. ...
    от billy_smoke в разделе Вопрошайка
  5. график функции в visual fortran
    нужно построить график функции в compaq visual fortran 6. нашла функцию plotp в imsl, почему-то ничего не строится (????), по оси x какие-то...
    от Anonymous в разделе Вопрошайка

Ваши права

  • Вы не можете создавать новые темы
  • Вы не можете отвечать в темах
  • Вы не можете прикреплять вложения
  • Вы не можете редактировать свои сообщения