Курс лекций по дисциплине Программирование на языке Фортран
 Скачать 0.98 Mb.
|
 ,  . При n=10, x=1, a=0.97 – D=0.1915155910-12.Интерполирование функций Интерполяционный многочлен Лагранжа:  или в более подробном виде: , где n – количество узлов интерполяции, xi – узлы интерполяции, f(xi) – значение функции в узле интерполяции. Если функция имеет сложный характер, то для вычисления ее значений для произвольного аргумента можно пользоваться более простыми интерполяционными формулами, построенными на значениях данной функции в некоторых узлах. Отметим, что интерполяционный многочлен будет иметь минимальную погрешность, если узлы не равноудалены друг от друга, а если они принимаются в виде: , m=1,...,n. ! File : Lagrang.for (l1.dat, l2.dat) ! Программа интерполирования с использованием интерполяционного ! многочлена Лагранжа ! n y(i)(x-x(0))*...*(x-x(i-1))*(x-x(i+1))*...*(x-x(n)) ! Ln(x)=СУМ ------------------------------------------------------------ ! i=0 (x(i)-x(0))*...*(x(i)-x(i-1))*(x(i)-x(i+1))*...*(x(i)-x(n)) ! n - количество узлов интерполирования; ! x(i)- узлы интерполяции; ! z - аргумент многочлена; ! F - значение многочлена при аргументе z real x(100),y(100),z real*8 F,L integer i,j,n character*12 fn write(5,*)' Введите имя файла исходных данных' read(5,'(a\)')fn open(2,file=fn) ! open(2,file='l2.dat') read(2,'(i20)')n ! write(5,*)n do i=1,n read(2,*)x(i),y(i) ! write(5,*)x(i),y(i) end do write(5,*)' Введите аргумент многочлена Лагранжа x' read(5,*)z F=0. do i=1,n L=1. do j=1,n if(i.ne.j) L=L*(z-x(j))/(x(i)-x(j)) end do ! Тест1 n=3 L=L*y(i) ! x 1 3 4 F=F+L ! y 12 4 6 L2(x)=2x¤-12x+22 end do ! При x=2, F=6 write(5,*)' F(',z,')=',F ! Тест2 n=4 read* ! x 0.41 1.55 2.67 3.84 end ! y 2.63 3.75 4.87 5.03 ! При x=1.91, F=4.153908736123 Формулы численного дифференцирования Формулы численного дифференцирования получены в результате дифференцирования интерполяционных формул:  ,  ,  . Более точные формулы: ,  .! File : differ.for ! Программа дифференцирования ф-ции f(x)=exp(2x) в точке по приближенным ! ф-лам f(x+h)-f(x-h) ! f'(x)=--------------- ! 2h ! ! f(x+h)-2f(x)+f(x-h) ! f''(x)=--------------------- ! h*h ! ! f(x-2h)-4f(x-h)+6f(x)-4f(x+h)+f(x+2h) ! f''''(x)=--------------------------------------- ! h**4 program Diff implicit none real x,h real*8 f,a,b write(*,*)'Введите шаг h' Read(*,*)h write(*,*)'Введите аргумент x ф-ции f(x)=exp(2x)' Read(*,*)x a=f(x-h) b=f(x+h) write(*,*)'Значение первой производной f''(x)=',(b-a)/(2.*h) write(*,*)' f''''(x)=',(b-2.*f(x)+a)/(h*h) read* ! Тест f'(x)=2exp(2x) end program ! При x=0 f'(0)=2 ! f''(x)=4exp(2x) ! При x=0 f''(0)=4 ! В программе при h=0.01, x=0 -> f'(x)=2.00013333599407, f''(x)=4.00013333510656 Function f(x) real x real*8 f f=exp(2.*x) end function ! Integral.f90 program Integral implicit none ! File: intgr3m.for ! Программа вычисления определенного интеграла 5 методами ! b ! Ї ! I(z)=¦f(x)dx, где f(x) - произвольная ф-ция ! ї ! a Integer n,j Integer*2 i1,i2,i3,i4 Real r Real*8 a,b Real*8 dx,i,k1,f,k3,xk,x,k2,f1,f2,f3,t1 Integer ng ! Кол-во узлов Гаусса 1<=ng<=3 Integer check ! Проверка размещения массива Real vg [ALLOCATABLE](:) ! Весовые коэффициенты Real ksig [ALLOCATABLE](:) ! Коэффициенты для определения узлов Гаусса в пределах интервала write(*,*)'Введите нижний и верхний пределы интегрирования a Read(*,*)a,b if (a.Gt.b) then k1=a ! Меняем местами пределы a=b b=k1 end if write(*,*)'Введите количество разбиений интервала (Num of devide interval)' write(*,*)'(для метода 4 это количество опытов)' Read(*,*)n dx=(b-a)/n ! Шаг разбиения ! 1. Формула трапеций I=0.0 k1=f(a) k3=0.5*dx xk=b-k3 do x=a,xk,dx k2=f(x+dx) I=I+k3*(k1+k2) k1=k2 end do write(*,*)'Значение интеграла по формуле трапеций (trapezoid) I=',I ! 2. Формула средних значений I=0.0 k1=a+0.5*dx ! f(x) вычисляется в n точках do x=k1,b,dx ! Погрешность eєM2*(b-a)**3/(24n¤) I=I+dx*f(x) ! M2=max¦f''(x)¦ при a end do write(*,*)'Значение интеграла по формуле средних (average) I=',I ! 3. Формула Симпсона I=0.0 k1=dx/6. ! f(x) вычисляется в 2n+1 точке k2=0.5*dx ! Погрешность eєM4*(b-a)**5/(180n**4) xk=b-k2 ! M4=max¦f''''(x)¦ при a f1=f(a) do x=a,xk,dx f2=f(x+k2) f3=f(x+dx) I=I+k1*(f1+4.*f2+f3) f1=f3 ! Значение справа на предыдущем шаге равно значению слева на последующем end do write(*,*)'Значение интеграла по формуле Симпсона (Sympson) I=',I ! 4. Метод Монте-Карло (применяется при многомерных интегралах) I=0.0 call gettim(i1,i2,i3,i4) call seed(i4) do j=1,n call random(r) x=a+(b-a)*r I=I+f(x) end do write(*,*)'Значение интеграла по методу Монте-Карло (Monte-Carlo) I=',I*(b-a)/n ! 5. Квадратурная формула Гаусса write(*,*)'Введите дополнительно количество узлов Гаусса (Num of Gausse nodes) (1-3)' Read(*,*)ng I=0.0 k2=0.5*dx allocate (vg(ng),ksig(ng),stat=check) ! Динамическое размещение массивов if(check.ne.0) stop 'Ошибка размещения массивов' Call CoeffGauss(ng,ksig,vg) ! Определение коэфф-тов Гаусса do j=1,ng ! Привидение коэффициента для [0,1] к интервалу [a,b] ksig(j)=ksig(j)*k2 end do k1=a+k2 ! Середина левого интервала do x=k1,b,dx ! Цикл по серединам интервалов do j=1,ng ! Цикл по узлам Гаусса k3=x+ksig(j) ! Координаты узлов Гаусса в пределах интервала I=I+f(k3)*vg(j) ! Значение функции*весовой коэффициент end do end do write(*,*)'Значение интеграла по формуле Гаусса (Gauss) I=',I*k2 write(*,*)'Количество узлов на интервале разбиения ng=',ng end program Integral Function f(x) Real*8 x,f f=x**4 return end subroutine CoeffGauss(ng,ksig,vg) ! Определение координат и весовых коэффициентов узлов Гаусса integer ng ! Кол-во узлов Гаусса 1<=ng<=3 integer check real vg(*) ! Весовые коэффициенты real ksig(*) ! Коэффициенты для определения узлов Гаусса в пределах интервала if (ng.lt.1.or.ng.gt.3) ng=3 select case (ng) case (1) ! При ng=1 формула Гаусса эквивалентна формуле средних ksig(1) = 0. vg(1) = 2. case (2) ksig(1) = -0.577350269189626 ! = -1# / Sqr(3#) = -0.57735 02691 89625 76451 ksig(2) = -ksig(1) vg(1) = 1. vg(2) = 1. case (3) ksig(1) = -0.774596669241483 ! = -Sqr(0.6) = -0.77459 66692 41483 37704 ksig(2) = 0. ksig(3) = -ksig(1) vg(1) = 0.555555555555556 ! = 5# / 9# = 0.(5) vg(2) = 0.888888888888889 ! = 8# / 9# = 0.(8) vg(3) = vg(1) end select end subroutine ! Тест: при f=x^4, a=0, b=1, n=4, ng=3 ! I1= 0.220703125000000 ! I2= 0.189697265625000 ! I3= 0.200032552083333 ! I4= случайное ! I5= 0.200000006729762 ДУ изгиба балки: , где y(x) – функция прогибов, M(x) – функция изгибающих моментов, EI – изгибная жесткость элемента. Прогиб изгибаемого элемента в сечении х можно также определить по формуле Мора: , где ! Deflection.f90 program Deflection implicit none ! Программа определения прогиба в заданном сечении x изгибаемого элемента при известной ! функции изгибающих моментов по формуле Мора с интегрированием по формуле Гаусса real*8 L ! Длина изгибаемого элемента, м integer n ! Количество участков разбиения integer ng ! Количество узлов Гаусса на участке разбиения real*8 x ! Сечение для определения прогиба real*8 dx ! Шаг разбиения real*8 EI ! Изгибная жесткость элемента real*8 M ! Изгибающий момент от фактической нагрузки real*8 M1 ! Изгибающий момент от единичной нагрузки real*8 q ! Внешняя равномерно-распределенная нагрузка real*8 f ! Прогиб real*8 t1,t2,Moment,OneForce real*8 xi [ALLOCATABLE](:) ! Координаты середин отрезков real*8 vg [ALLOCATABLE](:) ! Весовые коэффициенты real*8 ksig [ALLOCATABLE](:) ! Коэффициенты для определения узлов Гаусса в пределах интервала real*8 xg ! Координаты узлов Гаусса integer i,j L=6 q=4000 ! Н/м - внешняя нагрузка x=0.5*L ! Посередине !x=5.9 n=10 ng=3 EI=46.35E6 dx=L/n ! Шаг по длине элемента allocate (xi(n),vg(ng),ksig(ng)) ! Динамическое размещение массивов ! Вычисление xi(i)-координаты середины отрезка dx xi(1)=0.5*dx ! Координата первого сечения, м t1=xi(1) do i=2,n xi(i)=(i-1)*dx+t1 ! Определение координаты середины участка dx, м end do ! Вычисление упругого прогиба по формуле Мора (интегрирование по формуле Гаусса) ! (без учета влияния поперечной силы) f=0. call CoeffGauss(ng,ksig,vg) ! Определение коэфф-тов Гаусса do i=1,n t2=xi(i) do j=1,ng xg=t2+t1*ksig(j) ! Координаты узлов Гаусса по длине элемента M=Moment(q,L,xg) ! Вычисление момента в узле Гаусса M1=OneForce(xg,x,L) ! Изгибающий момент в узле Гаусса от единичной силы по направлению искомого прогиба f=f+M*M1*vg(j) ! Упругий прогиб от действия изгибающего момента end do end do f=1000.*f*t1/EI ! Прогиб от М, А b1=0.5*Deltax - вынесли постоянный множитель из цикла write(*,*)'Прогиб элемента в сечении х=',x,' равен ',f,' мм' end program Deflection ! Определение изгибающего момента в сечении х function Moment(q,L,x) real*8 Moment,q,L,x Moment=0.5*q*L*x-0.5*q*x*x ! Для шарнирно-опертой балки end function ! Определение изгибающего момента в сечении хi от единичной нагрузки в сечении х function OneForce(xi,x,L) real*8 OneForce,xi,x,Mx1,Qx1,L real*8 Ra ! Реакция левой опоры real*8 P1 ! Единичная сила Mx1=0. Qx1=0. P1=1. Ra=P1*(L-x)/L Qx1=Ra Mx1=Ra*xi if (xi.gt.x) then ! Входит сила Qx1=Qx1-P1 Mx1=Mx1-P1*(xi-x) end if OneForce=Mx1 end function subroutine CoeffGauss(ng,ksig,vg) integer ng ! Кол-во узлов Гаусса 1<=ng<=3 integer check real*8 vg(*) ! Весовые коэффициенты real*8 ksig(*) ! Коэффициенты для определения узлов Гаусса в пределах интервала if (ng.lt.1.or.ng.gt.3) ng=3 select case (ng) case (1) ! При ng=1 формула Гаусса эквивалентна формуле средних ksig(1) = 0. vg(1) = 2. case (2) ksig(1) = -0.577350269189626 ! = -1# / Sqr(3#) = -0.57735 02691 89625 76451 ksig(2) = -ksig(1) vg(1) = 1. vg(2) = 1. case (3) ksig(1) = -0.774596669241483 ! = -Sqr(0.6) = -0.77459 66692 41483 37704 ksig(2) = 0. ksig(3) = -ksig(1) vg(1) = 0.555555555555556 ! = 5# / 9# = 0.(5) vg(2) = 0.888888888888889 ! = 8# / 9# = 0.(8) vg(3) = vg(1) end select end subroutine ! Тест: однопролетная шарнирно-опертая балка длиной 6м с равномерно-распределенной нагрузкой ! q=4кН/м. Сечение прямоугольное 10х30см. Модуль упругости стали E=2.06*10^11Па ! Прогиб посередине: f=1.45631072344650 |