Проблема с нахождением корня нелинейного уравнения (Fortran 77) - Помощь студентам

reafantu · 02.03.2010, 17:40

Здравствуйте, в Fortran'e пока не силен и не могу никак найти свою ошибку. Программа находит корни нелинейного уравнения двумя способами: методом половинных отрезков и методом секущих, далее выводит график косинуса и моего.
Программа компилируется отлично, но после ввода данных она вылетает. Я нашел из-за чего, с помощью дебагера. Ошибка деления на ноль, но я не понимаю откуда она там и никак не могу исправить

Код всей программы:

Код:

PARAMETER (Na=100, Ni=10)
      Dimension XRAY(100),Y1RAY(100),Y2RAY(100)
      CHARACTER AXIS(2)*6,TITL(2)*22
      integer m,po,kar, ai
      real x(0:Ni), f(0:Ni)
      data (x(i), i = 1, Ni) /Ni*0/
      data (f(i), i = 1, Ni) /Ni*0/ 
      print'(''asin(x^m)-((1-x^j)^k'')'
      print'(''Vvedite* m,j,k: '',$) '
      read*,m,po,kar   
      print'(''Vvedite e,nmax: '',$) '
      read*,e,nmax
      a=0
      b=1
      n=0
      fa=FUNC(a,m,po,kar)
      
      do
      
      c=(a+b)/2
      fc=FUNC(c,m,po,kar)
      n=n+1
      s=a-b
      if ( s .lt. 0) s=s*(-1)
      if (E.gt. s) exit
      if (nmax.lt.n) exit
      if ((fa*fc) .gt. 0) then
      a=c
      fa=fc
      else
      b=c
      endif
      enddo
            
      imax=30
      i=1
      x(i-1)=0.2
      x(i)=x(0)+e*100
      ai=0
      f(i-1)=FUNC(x(i-1),m,po,kar)
      do
      if (ai .gt. imax) exit
      f(i)=FUNC(x(i),k,m,l)
      x(i+1)=x(i)-( (x(i)-x(i-1))/(f(i)-f(i-1)) )*(f(i))
      if (ABS(x(i+1)-x(i)) .lt. e) exit
      x(i-1)=x(i)
      f(i-1)=FUNC(x(i),k,m,l)
      x(i)=x(i+1)
      i=i+1
      ai=ai+1
      end do
      
      print*, 'chislo shagov=', ai
      print*, 'koren', x(i+1)
      print*, 'koren', x(i+1)
      print*, 'koren',(a+b)/2
      print*, 'n=', n
      read*
                 
      TITL(1)='FUNCTION'
      TITL(2)= 'SIN(X), COS(X)'     
      AXIS(1)='X-axis'
      AXIS(2)= 'Y-axis'      
      PI=3.1415926
      STEP=1./(Na-1)
      
      DO I=1,Na
      XRAY(I)=(I-1)*STEP
      Xa=XRAY(I)
      Y1RAY(I)=ASIN(Xa**m)-(1-xa**po)**kar
      Y2RAY(I)=COS(Xa)
      END DO  
         
      CALL GRAF_1(XRAY,Y1RAY,Y2RAY,Na,AXIS,TITL)
       
      WRITE(6,'(4X,A)') 'Press ENTER to finish ...'
      READ*   
      END
      
      function FUNC(x,m,po,kar)
      real m,po,kar
      FUNC=asin(x**m)-((1-x**po)**kar)
      return
      end


      SUBROUTINE GRAF_1(XRAY,Y1RAY,Y2RAY,Na,AXIS,TITL)
      CHARACTER AXIS(2)*6,TITL(2)*22
      Dimension XRAY(Na),Y1RAY(Na),Y2RAY(Na)                
      OY0=Y1RAY(1)
      OY1=Y1RAY(1)     
      OX0=XRAY(1)
      OX1=XRAY(1)        
      DO I=1,Na
      IF(OX0.GT.XRAY(I))     OX0=XRAY(I)
      IF(OX1.LT.XRAY(I))     OX1=XRAY(I)
      IF(OY0.GT.Y1RAY(I))    OY0=Y1RAY(I)
      IF(OY1.LT.Y1RAY(I))    OY1=Y1RAY(I)
      IF(OY0.GT.Y2RAY(I))    OY0=Y2RAY(I)
      IF(OY1.LT.Y2RAY(I))    OY1=Y2RAY(I)
      END DO
      OYH=(OY1-OY0)/4.
      OXH=(OX1-OX0)/4.
      CALL METAFL('XWIN')
      CALL SETPAG('DA4L') 
      CALL DISINI
      CALL PAGERA
      CALL HWFONT
      CALL AXSPOS(450,1800)
      CALL AXSLEN(2200,1200)
      CALL NAME(AXIS(1),'X')
      CALL NAME(AXIS(2),'Y')
      CALL LABDIG(1,'X')
      CALL TICKS(10,'XY')
      CALL TITLIN(TITL(1),1)
      CALL TITLIN(TITL(2),3)
      CALL GRAF(OX0,OX1,OX0,OXH,OY0, OY1,OY0,OYH)
      CALL TITLE
      CALL COLOR('RED')
      CALL CURVE(XRAY,Y1RAY,Na)
      CALL COLOR('GREEN')
      CALL CURVE(XRAY,Y2RAY,Na)
      CALL COLOR('FORE')
      CALL DASH
      CALL XAXGIT
      CALL DISFIN
       
      
      
      end

Ошибка именно в фрагменте:

Код:

do
      if (ai .gt. imax) exit
      f(i)=FUNC(x(i),k,m,l)
      x(i+1)=x(i)-( (x(i)-x(i-1))/(f(i)-f(i-1)) )*(f(i))
      if (ABS(x(i+1)-x(i)) .lt. e) exit
      x(i-1)=x(i)
      f(i-1)=FUNC(x(i),k,m,l)
      x(i)=x(i+1)
      i=i+1
      ai=ai+1
      end do

Serebro · 02.03.2010, 18:44

Цитата:

Сообщение от reafantu

Ошибка именно в фрагменте:

Код:

do
      if (ai .gt. imax) exit
      f(i)=FUNC(x(i),k,m,l)
      x(i+1)=x(i)-( (x(i)-x(i-1))/(f(i)-f(i-1)) )*(f(i))
      if (ABS(x(i+1)-x(i)) .lt. e) exit
      x(i-1)=x(i)
      f(i-1)=FUNC(x(i),k,m,l)
      x(i)=x(i+1)
      i=i+1
      ai=ai+1
      end do

Что такое l?

reafantu · 02.03.2010, 19:06

Неавнимательно делал, но проблема все равно не решилась

Код:

do
      if (ai .gt. imax) exit
      f(i)=FUNC(x(i),m,po,kar)
      x(i+1)=x(i)-( (x(i)-x(i-1))/(f(i)-f(i-1)) )*(f(i))
      if (ABS(x(i+1)-x(i)) .lt. e) exit
      x(i-1)=x(i)
      f(i-1)=FUNC(x(i),m,po,kar)
      x(i)=x(i+1)
      i=i+1
      ai=ai+1
      end do

Vago · 02.03.2010, 19:30

У Вас в основной программе m, po, kar (точнее, фактические параметры 2 - 4 при обращении к FUNC) - типа integer, а в функции они (точнее - формальные параметры 2-4) - real.

reafantu · 02.03.2010, 20:21

спасибо огромное, исправил. Но возникла новая ошибка...почему то не считает методом секущих корень....или считает но не точно (метод секущих это там где ошибка была)

Vago · 02.03.2010, 20:53

Цитата:

Сообщение от reafantu

....или считает но не точно

Чему m, po и kar равны?

Added 19:24 CET
1. Вы приближаетесь к корню слева и при этом у Вас a(0) - правее, чем a(1).

2. На результат это не влияет, но зачем Вы вообще храните последовательности x и f ?!

Код:

...
	imax=30
!	i=1
	xim1 = a ! 0.2
	xi = xim1 + E*100
	ai=0
	fim1 = FUNC( xim1, m, po, kar)
	do
	    if (ai .gt. imax) exit
	    fi = FUNC( xi , m, po, kar)
	    xip1 = xi - ( (xi-xim1) / (fi-fim1) )*(fi)
	    if (ABS(xip1-xi) .lt. e) exit
	    xim1 = xi
	    fim1 = fi   ! FUNC( xi, m, po, kar)
	    xi = xip1
!	    i=i+1
	    ai=ai+1
	end do
	print *, xi
...

reafantu · 03.03.2010, 16:30

m, po, kar по единице задавал. Спасибо за фрагмент программы, только почему то количество шагов всегда показывает единицу

P.S и заметил еще недоработку , по вашему фрагменту. Если брать эпсилон где-то 0.00001 тогда да , все считает, но если взять 0.1 то все...выкидывает

Vago · 03.03.2010, 19:42

Цитата:

Сообщение от reafantu

Спасибо за фрагмент программы

Пожалуйста.

Цитата:

Сообщение от reafantu

почему то количество шагов всегда показывает единицу

Вы перед методом секущих всегда гоняете метод дихотомии, а он у Вас портит a и b. Восстановите эти два значения перед началом работы по методу секущих, и всё встанет на свои места

Цитата:

Сообщение от reafantu

по вашему фрагменту

Вообще-то, это - Ваш фрагмент, и идея привязать x(1) к точности - Ваша идея

Я лишь подвинул x(0) на место и показал, как обойтись без массива.

Цитата:

Сообщение от reafantu

если взять 0.1 то все...выкидывает

Ну, естественно! X(1) оказался больше b. Можете попрактиковаться в Фортране, корректируя при необходимости x(1) перед тем как собственно приступать к поиску корня

03.03.2010, 16:30	#7
reafantu Регистрация: 02.03.2010 Сообщений: 4	m, po, kar по единице задавал. Спасибо за фрагмент программы, только почему то количество шагов всегда показывает единицу P.S и заметил еще недоработку , по вашему фрагменту. Если брать эпсилон где-то 0.00001 тогда да , все считает, но если взять 0.1 то все...выкидывает Последний раз редактировалось reafantu; 03.03.2010 в 16:37. Причина: добавление данных

Похожие темы
Тема	Автор	Раздел	Ответов	Последнее сообщение
Решение нелинейного уравнения и построение графика	xMass	Помощь студентам	5	29.10.2012 19:42
Численные методы. Решение нелинейного уравнения методом половинного деления.	gree	Помощь студентам	1	11.11.2009 18:36
Решение нелинейного уравнения методом Ньютона	Tina	Общие вопросы C/C++	2	04.06.2008 21:48

02.03.2010, 19:30	#4
Vago Форумчанин Регистрация: 15.01.2010 Сообщений: 948	У Вас в основной программе m, po, kar (точнее, фактические параметры 2 - 4 при обращении к FUNC) - типа integer, а в функции они (точнее - формальные параметры 2-4) - real.

02.03.2010, 20:21	#5
reafantu Регистрация: 02.03.2010 Сообщений: 4	спасибо огромное, исправил. Но возникла новая ошибка...почему то не считает методом секущих корень....или считает но не точно (метод секущих это там где ошибка была)