Здравствуйте уважаемые пользователи форума.
У меня возникла одна небольшая проблема.
Передо мной стояла следующая задача: найти зависимость скорости частицы(или любого другого объекта) в произвольной точке.
Значение в узлах сетки я считываю из предварительно данных мне файлов.
То есть у меня есть функция в трёх переменных в трёхмерном пространстве, которую нужно интерполировать с помощью полиномов Лагранжа
по заданным значениям в узлах сетки.
Случай интерполяции функции одной переменной не вызвал у меня затруднений. А вот здесь я немного сел в лужу.
http://stu.sernam.ru/book_dig_m.php?id=15
Это ссылка на формулу многочленов Лагранжа в двухмерном случае.
http://pmpu.ru/vf4/interpolation
А здесь случай произвольной размерности.
В трёхмерном случае узлы интерполяции, не должны лежать в одной плоскости. И здесь я поэтому запутался немного с циклами. Если кто-то найдёт ошибку или уже сталкивался с подобной проблей, то сообщите пожалуйста.
Уже месяц периодически мучаюсь с этой программой, которую писал сам, так что сразу извиняюсь за корявость кода.
Немного поясню суть метода Лагрнажа.
Пусть у нас есть область. В данном случае это прямоугольный параллелепипед. Мы разбиваем его на одинаковые маленькие прямоугольные паралеллепипеды. То есть получается такая трёхмерная сетка. В каждом узле сетки заданы значение функции. Теперь с помощью интерполяционных полиномов Лагранжа следует найти значение функции в произвольной точке.
Привожу код моей программы. Программа написана на языке Fortran.
Интерполяционные многочлены Лагранжа - Fortran
Код: Выделить всё
program time
call v_time
end
c------------------------------------------------------------------------
c Определение номера элемента, в который попадёт точка
subroutine place(xx,yy,zz,xstep,ystep,zstep)
real*8 zz(1), xx(1), yy(1), step, h, w, s
integer m, n, k, l, xstep(1), ystep(1), zstep(1)
c Задание точек, в которых нужно найти скорость жидкости
step = 0.1
c Размеры исходной области
h = 4 * acos(-1.0)/3
w = 1
s = 4 * acos(-1.0)
do m = 1, 1
xx(1) = 0.0
yy(1) = 0.0
zz(1) = 0.0
c xx(m + 1) = xx(m) + step
c yy(m + 1) = yy(m) + step
c zz(m + 1) = zz(m) + step
enddo
c Определение номера элемента по оси z, в который попадёт точка
do m = 1, 1
do n = 0, 24
if(zz(m).ge.(n*h/25) .and. (zz(m).le.((n + 1)*h/25)))
& then
zstep(m) = n + 1
endif
enddo
enddo
c Определение номера элемента по оси x, в который попадёт точка
do m = 1, 1
do l = 0, 29
if(xx(m).ge.(l*s/30) .and. (xx(m).le.((l + 1)*s/30)))
& then
xstep(m) = l + 1
endif
enddo
enddo
c Определение номера элемента по оси y, в который попадёт точка
do m = 1, 1
do k = 0, 9
if(yy(m).ge.(k*w/10) .and. (yy(m).le.((k + 1)*w/10)))
& then
ystep(m) = k + 1
endif
enddo
enddo
end
c--------------------------------------------------------------------------
c Вычисление базисных интерполяционных полиномов Лагранжа
subroutine lagrange(xx, yy, zz, x, y, z, l, i, j, k)
integer i, j, k, s, q, r
real*8 zz, xx, yy, x(8,8,8), y(8,8,8), z(8,8,8), l(8,8,8)
l = 1
do r = 1, 8
do s = 1, 8
do q = 1, 8
if(i.ne.r .and. j.ne.s .and. k.ne.q) then
l(i, j, k) = l(i, j, k) *
& * (xx - x(r,s,q))/(x(i,j,k) - x(r,s,q))
& * (yy - y(r,s,q))/(y(i,j,k) - y(r,s,q))
& * (zz - z(r,s,q))/(z(i,j,k) - z(r,s,q))
else
l(i, j, k) = l(i, j, k) * 1
endif
enddo
enddo
enddo
c write(*,*) l(i,j,k)
end
c--------------------------------------------------------------------------
c Сортировка точек в элементе с выделенной точкой
c для применения базисных Полиномов
subroutine polpoint(pox, poy, poz, num_step)
integer xstep(1), ystep(1), zstep(1)
integer l, n, num_step(512)
integer r, q, s
integer size_file, n1, el, gel, ex, ey, ez, i, j, k
real*8 x, y, z, zz(1), xx(1), yy(1), t(512)
real*8 u(512), v(512), pox(1,8,8,8)
real*8 poy(1,8,8,8), poz(1,8,8,8)
character(10) second_name
c Файл с расположением точек
second_name = 'Geom_total'
call place(xx,yy,zz,xstep,ystep,zstep)
open(unit = 1, file = second_name, form = 'unformatted',
& status='old')
read(1) size_file
c Я считываю данные из файла и записываю их так, чтобы применить метод Лагранжа.
do l = 1, size_file
read(1) n1, el, gel, ex, ey, ez, i, j, k, x,
& y, z
do n = 1, 1
if(ex.eq.xstep(n) .and. (ey-10).eq.ystep(n)) then
if(ez.eq.zstep(n)) then
do r = 1, 8
if(i.eq.r) then
do s = 1, 8
if(j.eq.s) then
do q = 1, 8
if(k.eq.q) then
t(q+8*(s-1)+64*(r-1)+512*(n-1)) = x
u(q+8*(s-1)+64*(r-1)+512*(n-1)) = y
v(q+8*(s-1)+64*(r-1)+512*(n-1)) = z
num_step(q+8*(s-1)+64*(r-1)+512*(n-1)) = l
endif
enddo
endif
enddo
endif
enddo
endif
endif
enddo
enddo
close(1)
c Отсортированные точки в спекиральном элементе
do n = 1, 1
do r = 1, 8
do s = 1, 8
do q = 1, 8
pox(n,r,s,q) = t(q+8*(s-1)+64*(r-1)+512*(n-1))
poy(n,r,s,q) = u(q+8*(s-1)+64*(r-1)+512*(n-1))
poz(n,r,s,q) = v(q+8*(s-1)+64*(r-1)+512*(n-1))
c write(3,*) pox(m,r,s,q), poy(m,r,s,q), poz(m,r,s,q)
enddo
enddo
enddo
enddo
end
c Сорртировка скоростей в узлах сетки, для применения базисных полиномов
Код: Выделить всё
subroutine velocity(num_step, third_name, v11, v22, v33)
character(16) third_name
real*8 vx, vy, vz, p
real*8 v1(512), v2(512), v3(512)
real*8 v11(1,8,8,8), v22(1,8,8,8), v33(1,8,8,8)
integer size_file, num_step(512), m, l, r, s, q
open(unit = 2, file = third_name, form='unformatted',
& status='old')
read(2) size_file
c Я считываю данные из файла и записываю их так, чтобы применить метод Лагранжа.
do l = 1, size_file
read(2) vx, vy, vz, p
do m = 1, 512
if(l.eq.num_step(m)) then
v1(m) = vx
v2(m) = vy
v3(m) = vz
endif
enddo
enddo
close(2)
c Отсортированный массив трёх компонет скорости
do m = 1, 1
do r = 1, 8
do s = 1, 8
do q = 1, 8
v11(m,r,s,q) = v1(q+8*(s-1)+64*(r-1)+512*(m-1))
v22(m,r,s,q) = v2(q+8*(s-1)+64*(r-1)+512*(m-1))
v33(m,r,s,q) = v3(q+8*(s-1)+64*(r-1)+512*(m-1))
c write(3,*) v11(m,r,s,q), v22(m,r,s,q), v33(m,r,s,q)
enddo
enddo
enddo
enddo
end
c--------------------------------------------------------------------------
c--------------------------------------------------------------------------
Вычисление среднего значения скорости для каждой точки
subroutine av_velocity(xx, yy, zz, pox, poy, poz, v11, v22, v33)
real*8 v11(1,8,8,8), v22(1,8,8,8), v33(1,8,8,8), lag(1,8,8,8)
integer i, j, k, n, xstep(1), ystep(1), zstep(1)
real*8 vx(1), vy(1), vz(1), xx(1), yy(1), zz(1)
real*8 pox(1,8,8,8), poy(1,8,8,8), poz(1,8,8,8)
call place(xx, yy, zz, xstep, ystep, zstep)
vx = 0
vy = 0
vz = 0
lag = 1
do n = 1, 1
do i = 1, 8
do j = 1, 8
do k = 1, 8
c write(*,*) v11(n,i,j,k), v22(n,i,j,k), v33(n,i,j,k)
c Применние интерполяционных полиномов Лагранжа
call lagrange(xx(n),yy(n),zz(n),pox(n,:,:, :) ,poy(n,:,:, :) ,
& poz(n,:,:, :) , lag(n,:,:, :) , i, j, k)
vx(n) = vx(n) + v11(n,i,j,k) * lag(n,i,j,k)
vy(n) = vy(n) + v22(n,i,j,k) * lag(n,i,j,k)
vz(n) = vz(n) + v33(n,i,j,k) * lag(n,i,j,k)
enddo
enddo
enddo
enddo
call place(xx, yy, zz, xstep, ystep, zstep)
do n = 1, 1
c write(3,*) xx(n), yy(n), zz(n), vx(n), vy(n), vz(n)
enddo
end
c----------------------------------------------------------------------------
c Нахождение зависимости эволюции средних скоростей
c частиц в зависимости отвремени
subroutine v_time
integer i, num_step(512)
character(6) first_name
character(16) third_name
real*8 v11(1,8,8,8), v22(1,8,8,8), v33(1,8,8,8), poz(1,8,8,8)
real*8 xx(1), yy(1), zz(1), pox(1,8,8,8), poy(1,8,8,8)
c call polpoint(lag, num_step)
third_name = 'Velo_total_00001'
first_name = '5_file'
open(unit = 3, file = first_name)
call polpoint(pox, poy, poz, num_step)
do i = 1, 1
write(third_name(12:16),'(I5.5)') i
call velocity(num_step, third_name, v11, v22, v33)
call av_velocity(xx, yy, zz, pox, poy, poz, v11, v22, v33)
write(3,*) " "
enddo
close(3)
end
Зависимость от чего?найти зависимость скорости частицы(или любого другого объекта) в произвольной точке.
Немного описана суть интерполяции, про метод Лагранжа не сказано ни слова. Кроме того, скорость - это первая производная от пути. То есть не исключено применение интерполяционных полиномов Эрмита.Немного поясню суть метода Лагрнажа.
Кода нет.Привожу код моей программы.
Даже самый дурацкий замысел можно воплотить мастерски
Спасибо за уточнения.
1) Нужно найти скорость частицы в зависимости от времени.
2) Тут я неправильно сформулировал, это суть интерполяции. Про метод Лагранжа можно почитать по ссылкам выше.
3) У вас такие правила по ограничению количества символов, поэтому я и подгружал по частям.
1) Нужно найти скорость частицы в зависимости от времени.
2) Тут я неправильно сформулировал, это суть интерполяции. Про метод Лагранжа можно почитать по ссылкам выше.
3) У вас такие правила по ограничению количества символов, поэтому я и подгружал по частям.