newton para sistemas não lineares

Esta é uma pré-visualização de arquivo. Entre para ver o arquivo original

program newton
implicit none
real::E,j1,j2,n1,n2,n3,n4,n5,n6,m,z1,z2
real, dimension(100)::x,y,v1,v2,j3,j4,f1,f2,t1,t2
real, dimension(2)::s2,s1
integer::i
E=0.000000000001
x(1)=1.0
y(1)=5.0
!Passo 1: Calcular F(x(i)) e J(x(i));
i=1
20 f1=x(i)+y(i)-3.0
 f2=x(i)**2+y(i)**2-9.0
 J1=1.0
 j2=1.0
 j3=2*x(i)
 j4=2*y(i)
 t1(i)=m(f1(i),f2(i))
if (t1(i)<E) then
write(*,*)i,t1(i),x(i),y(i)
else
write(*,*)
write(*,*)
!Passo 3: obtenho V(i), solução do sistema linear: J(x(i))*V(i)=F(x(i));
 s1(1)=0.0
 s1(2)=0.0
 n1=f1(i)
 n2=j2
 n3=j1
 n4=f2(i)
 n5=j3(i)
 n6=j4(i)
 
!write(*,*)n1,n2,n3
!write(*,*)n4,n5,n6
 call jacobi(s1,s2,n1,n2,n3,n4,n5,n6)
 
v1(i)=s2(1)
v2(i)=s2(2)
write(*,*)i,v1(i),v2(i)
 x(i+1)=x(i)+v1(i)
 y(i+1)=y(i)+v2(i)
z1=x(i+1)-x(i)
z2=y(i+1)-y(i)
t2(i)=m(z1,z2)
if (t2(i)<E) then
write(*,*)
write(*,*)i,t2(i),x(i+1),y(i+1)
else
i=i+1
go to 20
continue
end if
end if
!Volte ao passo 1.
end program newton
! Subrotina para resolver o sistema linear J(x(i))*v(i)=F(x(i)), onde j(x)= jacobiana b= f(i)
subroutine jacobi(s1,s2,n1,n2,n3,n4,n5,n6)
 implicit none
real,intent(in)::n1,n2,n3,n4,n5,n6
!s1,n1,n2,n3,n4,n5,n6 são parâmetros de entrada e s2 é um parâmetro de saída da subrotina. s1 e s2 são vetores de R² e os outros são reais
real, dimension(2),intent(in)::s1
real, dimension(2), intent(out)::s2
 integer :: i
 real, dimension(100) :: x,y
x(1) = s1(1)
y(1) = s1(2)
do i = 1, 99
 x(i+1) = (-n1-n2*y(i))/n3
 y(i+1) = (-n4-n5*x(i))/n6
 
end do
s2(1)=x(100)
s2(2)=y(100)
return
end subroutine jacobi
function m(x,y)
real::x,y,m
m=max(abs(x),abs(y))
return
end function m