Usuário:Manoel
De WikiLICC
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