Usuário:Manoel

Página Wiki do Manoel.

Testar quando n=3?

Olhar (inv(Q)*S' )?

Codigo convertido

program trabalho1_1       !function [u]=trabalho1_1(n,xp,yp)
 implicit none
!clear
!n=4;
!a1=1;
!a2=2;

!tic;

 integer,parameter       ::  n=4
 integer                 ::  i,j
 real*8                  ::  xp,yp,fn

 real*8,dimension(n)     ::  x,y
 real*8,dimension(n,n)   ::  XX,YY,deltaX,deltaY,FV
 real*8,dimension(n*n,1) ::  Sx_1,Sy_2
 real*8,dimension(n*n,1) ::  F,FVV,deltaF
 real*8,dimension(n*n,9) ::  S
 real*8,dimension(9,n*n) ::  SdF,ST
 real*8,dimension(9,9)   ::  Q,invQ
 real*8,dimension(9,1)   ::  f_nx1
    
 xp=0.2
 yp=0.3


 !x=linspace(0,1,n);       !partiçao do espaço com n pontos
 do i=1,n
   x(i)=real(i)/n
 end do
  
 XX=0.D0                   !X=[zeros(n,n)]

 do i=1,n                  !for i=1:n
    XX(:,i)=x(i)           !X(:,i)=x(1,i)        ! as n's linhas igual a primeira 
 end do
 !pause
 !=================================================================
 do i=1,n                  !y=linspace(0,1,n);          
   y(i)=real(i)/n          !partiçao do espaço com n pontos
 end do

 YY=0.D0                   ! Y=[zeros(n,n)] !só com zeros dá
 do i=1,n                  !  for i=1:n
    YY(i,1)=y(i)           !  Y(i,:)=y(1,i)        ! as n's colunas igual a primeira
 end do
 !pause
 !=================================================================
 FN=xp**2+yp**2                       !função nodo onde a1 e a2 são coordenada dos nodos
 do i=1,n                             !for i=1:n
    do j=1,n                          !  for j=1:n
       FV(i,j)=XX(i,j)**2+YY(i,j)**2  !Função da vizinhança 
    end do
 end do

 FVV=reshape(FV, (/n*n,1/))           !FVV=FV(:)    ! todos os nos vizinhos em forma de coluna

 F     =0.D0                          !F=zeros(n*n,1)            !F tem n*n linhas e 1 coluna  
 F     =FN                            !F(:,1)=FN              !toda primeira coluna é igual a Função nodo
 DeltaF=FVV-F                         !coluna Delta 

 !=================================================================
 DeltaX=XX-xp                         !DeltaX=X-xp        ! todos oa itens da matriz menos nodo central
 Sx_1=reshape(DeltaX, (/n*n,1/))      !Sx_1=DeltaX(:)      ! matriz virou em coluna (vetor)
 !=================================================================
 DeltaY=YY-yp                         ! todos oa itens da matriz menos nodo central
 Sy_2=reshape(DeltaY, (/n*n,1/))      !Sy_2=DeltaY(:)      ! matriz virou em coluna (vetor)
 !=================================================================
 S(:,1)=Sx_1(:,1)
 S(:,2)=Sy_2(:,1)
 do i=1,n*n                              !for i=1:n*n;
    S(i,3)=(Sx_1(i,1)**2)/2.D0            ! S(i,3)=(Sx_1(i,1)^2)/2
    S(i,4)= Sx_1(i,1)*Sy_2(i,1)
    S(i,5)=(Sy_2(i,1)**2)/2.D0
    S(i,6)=(Sx_1(i,1)**3)/6.D0
    S(i,7)=((Sx_1(i,1)**2) *Sy_2(i,1))/2.D0
    S(i,8)=(Sx_1(i,1)*Sy_2(i,1)**2)/2.D0
    S(i,9)=(Sy_2(i,1)**3)/6.D0
 end do

 !S(:,3)= (Sx_1(:,1).^2)/2
 !
 !=================================================================
 ST= transpose(S)
 Q = matmul(ST,S)       !Q=S'*S           ! 1 multiplicar pela transposta para achar matriz quadrada
 !NAO TEM                cond(Q)
 !W1=inv(Q)

 !invQ = inv( Q )
 !FALTA CALCULAR A INVERSA DE Q
 invQ = Q !Apenas para continuar a execucao
  
 SdF  = matmul( ST,DeltaF )
 f_nx1= matmul(invQ ,SdF  )     !f_nx1=inv(Q)*(S'*DeltaF)
 print *,f_nx1
end program
!whos            !mostra os vetores
!toc