Форум программистов [Powered by Invision Power Board]

Версия для печати темы
Нажмите сюда для просмотра этой темы в оригинальном формате

Форум программистов > Fortran > [General] Метод бисекции

Автор: FFFU 18.11.2010, 00:33

Здравствуйте!
Нам в институте дали задание написать программу для нахождения корней функции методом деления пополам.

Код


program func
dimension x(101), v(101), w(101), u(101)
real x, y, t, s
pi=3.141592653589793
a=0; b=1; n=100
c=2; l=1; m=3; t=5; q=2; s=0.5
h=(b-a)/n

do i=1, n+1
 x(i)=a+h*(i-1)
 v(i)=c*sin(l*(pi/2)+pi*(x(i)**m))
 w(i)=-t*(x(i)**q)+s
 u(i)=v(i)+w(i)
 print*, i, x(i), u(i)
enddo

end

надо найти точку пересечения функции u(i) с осью абсцисс, то есть ноль функции.
Я почитал в википедии в чем суть этого метода и как его нужно примерно делать. Но я не понимаю, как это сделать с функцией u(i). По-моему, она как бы зависит он счетчика i. И как тогда задать начало и конец? u(1) и u(101)? Но ведь это всего два числа, как задать, что это именно отрезок? Как разделить его пополам? Я бы мог понять, если бы нужно было сделать один большой массив и найти в нем число, в котором функция ближе всего к нулю, но тут требуется вроде как что-то другое, а что я понять не могу.
Помогите, пожалуйста.

Автор: FCM 18.11.2010, 13:08

Функцию надо вынести в отдельную фортран-функцию.
Упрощенно будет выглядеть примерно так.

Код


program zero
     implicit none
     interface
         real function f(x)
            real, intent(in) :: x
         end function f
     end interface
     real :: a, b, x, eps = 1.e-6
11 write(*,*) 'enter two points ( f(a)*f(b) must be <= 0 )'
     read(*,*) a, b
     if( f(a)*f(b) > 0 ) go to 11
    
   ... здесь алгоритм нахождения 0 методом бисекции

     write(*,*) 'x = ', x, ';   f(x) = ', f(x)
end

real function f(x)
   implicit none
   real, intent(in) :: x
   real, parameter :: pi = 3.141592653589793
   real :: c, l, m, t, q, s, v, w
   c=2; l=1; m=3; t=5; q=2; s=0.5
   v = c*sin( l*(pi/2) + pi*(x**m) )
   w = -t*(x**q) + s
   f = v + w
end

PS/ Метод решения также можно вынести в отдельную фортран-функцию

Автор: FFFU 24.11.2010, 01:06

Простите, что не отвечал так долго, готовился к контрольной по химии.
Вот алгоритм, который я честно взял из http://ru.wikipedia.org/wiki/%D0%9C%D0%B5%D1%82%D0%BE%D0%B4_%D0%B1%D0%B8%D1%81%D0%B5%D0%BA%D1%86%D0%B8%D0%B8#.D0.9F.D1.80.D0.BE.D0.B3.D1.80.D0.B0.D0.BC.D0.BC.D0.BD.D1.8B.D0.B9_.D0.BA.D0.BE.D0.B4 smile

Пока без условия про крайние точки.

Код


      program zero
      implicit none
      interface
        real function f(x)
         real, intent(in) :: x
        end function f
      end interface
      real :: a,b,x,r,eps=1.e-6
   11 write(*,*) ' Введите a и b. f(a)*f(b) < 0'
       read (*,*) a,b
       if (f(a)*f(b).gt.0) go to 11
        do while (f(r).ne.0.and.ABS(b-a).gt.eps)
         r=(a+b)/2
         if(f(a)*f(r).lt.0) b=r
         if(f(r)*f(b).lt.0) a=r
        enddo
       write(*,*) 'x= ',x, 'f(x)= ', f(x)
      end

В программе, которая была в первом сообщении было видно, что 0 находится где-то между 0.62 и 0.63 икса.
А нынешняя выдает

Код



 Введите a и b. f(a)*f(b) < 0
0,1
 x=   1.4012985E-45 f(x)=    2.500000

Первое число слишком маленькое, а второе это максимум функции на отрезке [0;1].
Подскажите, пожалуйста, где я ошибся.

Автор: FCM 24.11.2010, 07:21

Приведу свой вариант

Код


program func
   implicit none
   interface
      real function f(x)
         real, intent(in) :: x
      end function f
   end interface
   real :: a, b, x, eps = 1.e-6
11 write(*,*) 'enter two points ( f(a)*f(b) must be <= 0 )'
   read(*,*) a, b
   if( f(a)*f(b) > 0 ) go to 11
   do while( abs(a-b) > eps )
       x = (a+b)/2; 
       if( f(a)*f(x) > 0 )    then
         a = x
       else
         b = x
       end if
   end do
   write(*,*) 'x = ', x, ';   f(x) = ', f(x)
end

real function f(x)
   implicit none
   real, intent(in) :: x
   real, parameter :: pi = 3.141592653589793
   real :: c, l, m, t, q, s, v, w
   c=2; l=1; m=3; t=5; q=2; s=0.5
   v = c*sin( l*(pi/2) + pi*(x**m) )
   w = -t*(x**q) + s
   f = v + w
end

Результат на отрезке [0,1]
x = 0.62384319 ; f(x) = -9.77516174E-06

Автор: FFFU 25.11.2010, 14:28

Большое спасибо за помощь!
Все заработало. Ошибка оказалась в строке
write(*,*) 'x= ',x, 'f(x)= ', f(x)
у меня была переменная r, а я ее забыл подставить