Resolución de Sistemas de Ecuaciones Lineales y Cálculo Numérico

Problema 41

if(A.dim1()==0 || A.dim1()!=A.dim2()) return 0.;

if(A.dim1()==1) return A[0][0];

real determinante=0.;
  for(int k=0;k    Array2D B(A.dim1()-1,A.dim1()-1);
    for(int i=0;i      for(int j=0;j        if(j        else  B[i][j]=A[i+1][j+1];
      }
    }
    if(k%2==0) determinante+=A[0][k]*mn_determinante_recursivo(B);
    else determinante-=A[0][k]*mn_determinante_recursivo(B);
  }
  return determinante;

Problema 42

Array2D A=A_original.copy();
  Array1D b=b_original.copy();

  if(A.dim1()!=A.dim2() || A.dim1()!=b.dim() || b.dim()==0) return Array1D();

  for(int k=0;k    int kmax=max_pos(A,k);
    if(A[kmax][k]==0.) return Array1D();
    if(kmax!=k){
      for(int j=k;j        mn_pivotar(A[k][j],A[kmax][j]);
      }
      mn_pivotar(b[kmax],b[k]);
    }
    for(int j=k+1;j          real mul=-A[j][k]/A[k][k];
          A[j][k]=0.;
          for(int n=k+1;n          b[j]+=mul*b[k];
    }
  }

  return mn_remonte(A,b);

Problema 43

Array2D A=A_original.copy();
  Array1D b=b_original.copy();

  if(A.dim1()!=A.dim2() || A.dim1()!=b.dim() || b.dim()==0) return Array1D();

  Array1D  piv(b.dim());
  for(int k=0;k

  for(int k=0;k    int kmax=max_pos(A,k,piv);
    if(A[piv[kmax]][k]==0.) return Array1D();
    if(kmax!=k) mn_pivotar(piv[k],piv[kmax]);

    for(int j=k+1;j          real mul=-A[piv[j]][k]/A[piv[k]][k];
          A[piv[j]][k]=0.;
          for(int n=k+1;n          b[piv[j]]+=mul*b[piv[k]];
    }
  }

  return mn_remonte(A,b,piv);

Problema 44

Array2D A=A_original.copy();
  if(A.dim1()!=A.dim2() || A.dim1()==0) return Array2D();

  Array2D B(A.dim1(),A.dim2(),0.);
  for(int k=0;k

  for(int k=0;k    int kmax=max_pos(A,k);
    if(A[kmax][k]==0.) return Array2D();

    if(kmax!=k){
      for(int i=0;i        mn_pivotar(A[k][i],A[kmax][i]);
        mn_pivotar(B[k][i],B[kmax][i]);
      }
    }

    for(int j=k+1;j          real mul=-A[j][k]/A[k][k];
          A[j][k]=0.;
          for(int n=k+1;n          for(int n=0;n    }

Problema 45

  if(M.dim1()==0 || M.dim1()!=M.dim2()) return 0.;

  Array2D A=M.copy();
  real determinante;

  int Npiv=0;
  for(int k=0;k    int kmax=max_pos(A,k);
    if(A[kmax][k]==0.) return 0.;

    if(kmax!=k){
      Npiv++;
      for(int i=0;i    }

    for(int j=k+1;j        real mul=-A[j][k]/A[k][k];
        A[j][k]=0.;
        for(int n=k+1;n    }
  }

  determinante=1;
  for(int k=0;k

  if(Npiv%2==1) return(-determinante);
  return determinante;

Problema 46

1)

   if(A.dim1()!=A.dim2() ) return( Array2D());
   int N=A.dim1();

   Array2D B(N,N);

   for (int i=0;i      real Sum = 0.;
      for (int k=0;k      if (A[i][i]());
      B[i][i]=sqrt(A[i][i]-Sum);
      for  (int j=i+1;j      Sum = 0;
         if (B[i][i]==0) return( Array2D());
         for (int k=0;k         B[j][i]=(A[j][i]-Sum)/B[i][i];
         B[i][j]= B[j][i];
      }
   }
   return( B );

2)

int N=A.dim1();
    Array2D CH = mn_cholesky_factorization (A);
    if (CH.dim1()==0) return (Array1D() );
    Array1D z=mn_descenso (CH,b);
    if (z.dim()==0) return (Array1D() );
    Array1D u = mn_remonte (CH,z);
    if (u.dim()==0) return (Array1D() );
    return(u);

Problema 47

1)

   int N=A.dim1();
   L=Array2D(N,N,0.);
   U=Array2D(N,N,0.);

   if(N!=A.dim2() || N!=L.dim1() || N!=L.dim2() || N!=U.dim1() || N!=U.dim2()) return(-1);

   for(int i=0;i      L[i][i]=1;
      real Sum = 0.;
      for(int k=0;k      U[i][i]=A[i][i]-Sum;
      if(U[i][i]==0.) return(-1);
      for(int j=i+1;j       Sum = 0;
        for (int k=0;k        U[i][j]=A[i][j]-Sum;
        Sum = 0;
        for (int k=0;k        L[j][i]=(A[j][i]-Sum)/U[i][i];
      }
   }
   return(0);

2)

Array1D c=mn_descenso(L,b);
return mn_remonte(U,c);

Problema 48

1)

  if(l.dim()!=a.dim() || l.dim()!=(m.dim()+1) || u.dim()!=m.dim() || u.dim()!=b.dim() || u.dim()!=c.dim()) return(-1);

  l[0]=a[0];
  if(l[0]==0) return(-2);
  u[0]=b[0]/l[0];
  for(int i=1;i    m[i-1]=c[i-1];
    l[i]=a[i]-m[i-1]*u[i-1];
    if(l[i]==0) return(-2);
    u[i]=b[i]/l[i];
  }
  m[m.dim()-1]=c[c.dim()-1];
  l[a.dim()-1]=a[a.dim()-1]-m[a.dim()-2]*u[a.dim()-2];
  return(0);

2)

  Array1D l(a.dim());
  Array1D m(c.dim());
  Array1D u(b.dim());

  int error=crout_descomposicion(a,b,c,l,m,u);
  if(error());

  Array1D z=crout_descenso (l,m,t);
  if(z.dim()==0) return(Array1D());

  return crout_remonte (u,z);

Problema 61

  Array1D autovalores(N);
  Autovectores=Array2D(N,N,0.);
  for(int i=0;i

  for (Niter=0;Niter    int p=0;
    int q=1;
    real Max = 0;
    for (int i=0;i      for (int j=i+1;j        if (fabs(A[i][j])>Max){
          Max=fabs(A[i][j]);
          p=i;
          q=j;
        }
      }
    }
    if (Max      for(int l=0;l      return (autovalores);
    }

    real alfa=0.5*atan( (2.*A[p][q])/(A[q][q]-A[p][p]) );
    real co=cos(alfa);
    real si=sin(alfa);

    for(int k=0;k      real temp= Autovectores[k][p];
      Autovectores[k][p]=co*temp-si*Autovectores[k][q];
      Autovectores[k][q]=co*Autovectores[k][q]+si*temp;
      if(k==p || k==q) continue;
      temp=A[p][k];
      A[k][p]=A[p][k]=co*temp-si*A[q][k];
      A[k][q]=A[q][k]=co*A[q][k]+si*temp;
    }
    real App=A[p][p];
    A[p][p]=co*(co*App-si*A[p][q])-si*(co*A[p][q]-si*A[q][q]);
    A[q][q]=co*(co*A[q][q]+si*A[p][q])+si*(co*A[p][q]+si*App);
    A[p][q]=A[q][p]=0.;

  }

  return(Array1D());

Problema 62

  l_max=mn_norma_euclidea(u_max);
  u_max=u_max/l_max;
  for(int i=0;i    Array1D u_new=A*u_max;
    real l_new=mn_norma_euclidea(u_new);
    if( l_new==0) return(-2);
    if(mn_producto_escalar(u_new,u_max)    else l_max=l_new;
    u_new=u_new/l_max;
    if( mn_norma_euclidea(u_new-u_max)      u_new=u_max;
      return(i);
    }
    u_max=u_new;
  }
  return(-3);

Problema 63

  l_min=mn_norma_euclidea(u_min);
  u_min=u_min/l_min;
  for(int i=0;i    Array1D u_new=mn_gauss(A,u_min);
    if(u_new.dim()==0) return -1;
    real l_new=mn_norma_euclidea(u_new);
    if( l_new==0) return(-2);
    if(mn_producto_escalar(u_new,u_min)    else l_min=l_new;
    u_new=u_new/l_min;
    if( mn_norma_euclidea(u_new-u_min)      u_new=u_min;
      l_min=1./l_min;
      return(i);
    }
    u_min=u_new;
  }
  return(-3);

  for(int k=0;k    A[k][k]-=l_min;
  }
  real l_min0=l_min;
  l_min=u_min.norma();
  u_min=u_min/l_min;
  for(int i=0;i    Array1D u_new=mn_gauss(A,u_min);
    if(u_new.dim()==0) return -1;
    real l_new=u_new.norma();
    if( l_new==0) return(-2);
    if(u_new*u_min    else l_min=l_new;
    u_new=u_new/l_min;
    if((u_new-u_min).norma()      u_new=u_min;
      l_min=l_min0+1./l_min;
      return(i);
    }
    u_min=u_new;
  }
  return(-3);

Problema 65

  real lambda=1.;

  for(int k=0;k    real Fu=F(u);
    Array1D Gu=mn_gradiente(F,u,h);

    Array1D v=u-Gu*lambda;
    real Fv=F(v);

    while(Fv>=Fu){
      if(mn_average_distance(u,v)      lambda/=10.;
      v=u-Gu*lambda;
      Fv=F(v);
    }

    if( mn_distancia(Fu,Fv)      u=v;
      return k;
    }

    lambda*=10;
    u=v;
    Fu=Fv;
  }

  return Nmax;

Problema 66

  Array1D fu=f(u);
  for(int i=0;i    if(fu.norm()    Array2D Derivada=matriz_derivada(f,u);
    Array1D sol=mn_gauss(Derivada,fu*(-1.));
    if(sol.dim()==0) return(-1);
    if(sol.norm()    u=u+sol;
    fu=f(u);

  }
  return(-1);

Problema 67

real lambda=1.;

  real Error=Error_p(u,x,y,p);
  Array1D fu=Grad_Error_p(u,x,y,p,h);
  if(fu.dim()==0) return(-1);

  for(int i=0;i    if(fu.norm()==0. || lambda    Array2D Derivada=matriz_derivada(u,x,y,p,h)*lambda;
    if(Derivada.dim1()==0) return(i);
    for(int k=0;k    Array1D sol=mn_gauss(Derivada,fu*(-lambda));
    if(sol.dim()==0){
      lambda/=10.;
      continue;
    }
    if(sol.norm()    Array1D u_new=u+sol;
    real ErrorNew=Error_p(u_new,x,y,p);
    if(ErrorNew      fu=Grad_Error_p(u_new,x,y,p,h);
      if(fu.dim()==0) return(i);
      u=u_new;
      Error=ErrorNew;
      lambda*=10.;
    }
    else{
      lambda/=10.;
    }
  }
  return(Nmax);

Problema 68

if(A.dim1()!=A.dim2() || A.dim1()!=b.dim() || u.dim()!=b.dim() || b.dim()==0){
    return -1;
  }
  for(int k=0;k    if(A[k][k]==0.) return -1;
  }

  for(int n=1;n    real error=0.;
    if(n%2==0){
      for (int i=0;i        real suma=0;
        for(int j=0;j          if(i!=j) suma = suma - A[i][j]*u[j];
        }
        suma = (suma + b[i])/A[i][i];
        error+=mn_distancia(suma,u[i]);
        u[i] = suma;
      }
    }
    else{
      for (int i=b.dim()-1;i>=0;i--){
        real suma=0;
        for(int j=0;j          if(i!=j) suma = suma - A[i][j]*u[j];
        }
        suma = (suma + b[i])/A[i][i];
        error+=mn_distancia(suma,u[i]);
        u[i] = suma;
      }
    }
    if (error/b.dim()  }
  return -1;

Problema 69

if(A.dim1()!=A.dim2() || A.dim1()!=b.dim() || u.dim()!=b.dim() || b.dim()==0){
     return -1;
    }
    for(int k=0;k      if(A[k][k]==0) return -1;
    }
    for(int n=1;n       real error=0;
       for (int i=0;i           real suma=0;
           for(int j=0;j              if(i!=j){
                 suma = suma - A[i][j]*u[j];
              }
           }
           suma = w*(suma + b[i])/A[i][i]+(1-w)*u[i];
           error+=mn_distancia(suma,u[i]);
           u[i] = suma;
       }
       if (error/b.dim()    }
    return -1;

Problema 610

    if(A.dim()!=b.dim() || u.dim()!=b.dim() ){
     printf("relajacion() : arrays of different dimensions \n");
    }
    int dim=b.dim1();

    int i,j,n;
    real temp2;
    for(n=1;n    {
       Array1D u2=u.copy();
       for (i=0;i           temp2=0;
           for(j=1;j              temp2 = temp2 - A[i][j]*u[J[i][j]];
           }
           temp2 = temp2 + b[i];
           u[i] = w*temp2/A[i][0]+(1-w)*u[i];
       }
       if (mn_error_vectores(u,u2)    }
    return -1;

Problema 611

1)

  real norma=0;
  for(int j=0;j    real temp=0;
    for(int i=0;i      temp+=fabs(A[i][j]);
    }
    if(temp>norma) norma=temp;
  }
  return norma;

2)

  real norma=0;
  for(int i=0;i    real temp=0;
    for(int j=0;j      temp+=fabs(A[i][j]);
    }
    if(temp>norma) norma=temp;
  }
  return norma;

3)

  Array2D B=A.transpose()*A;
  return sqrt(mn_autovalor_maximo(B));

Problema 35

  Array1D h(x.dim()-1);
  Array1D M(x.dim(),1.);
  Array1D B(x.dim());
  Array1D U(x.dim()-1,0.);
  Array1D L(x.dim()-1,0.);

  for(int k=0;k

  B[0]=c0;
  B[B.dim()-1]=cN;

  for(int k=1;k    L[k-1]=h[k-1];
    M[k]=2*(h[k-1]+h[k]);
    U[k]=h[k];
    B[k]=3.*(f[k+1]-f[k])/h[k] - 3.*(f[k]-f[k-1])/h[k-1];
  }
  c=solucion_sistema(L,M,U,B);

  a=Array1D(x.dim()-1);
  b=Array1D(x.dim()-1);
  d=Array1D(x.dim()-1);

  for(int k=0;k    a[k]=f[k];
    d[k]=(c[k+1]-c[k])/(3*h[k]);
    b[k]=(f[k+1]-f[k])/h[k]-h[k]*(2*c[k]+c[k+1])/3.;
  }

  return 0;
}

Problema 27

1)

  if(a.dim()==1) return Array1D();
  Array1D b(a.dim()-1);
  for(int k=0;k    b[k]=(k+1.)*a[k+1];
  }
  return b;

2)

  if(a[a.dim()-1]==0.) return d;

  real max_c=mn_abs(a[0]);
  for(int k=1;k    if(mn_abs(a[k])>max_c) max_c=mn_abs(a[k]);
  }
  real A=1+max_c/mn_abs(a[a.dim()-1]);

  Array1D ceros;

  Array1D I(1,-A);
  for(int k=0;k  I.push_back(A);

  for(int k=0;k    if(evaluar_polinomio(a,I[k])*evaluar_polinomio(a,I[k+1])      ceros.push_back(calculo_cero_en_intervalo(a,I[k],I[k+1],TOL));
    }
  }

  return ceros;

3)

  Array1D ceros;
  for(int k=a.dim()-2;k>=0;k–){
    Array1D d=a.copy();
    for(int n=0;n    ceros=ceros_polinomio_desde_ceros_derivada(d,ceros,TOL);
  }
  return ceros;

51)

   if(F.dim1()();
  Array2D M(F.dim1(),F.dim2(),0.);
  for(int i=1;i    for(int j=1;j      for(int k=-1;k        for(int l=-1;l          M[i][j]+=m[k+1][l+1]*F[k+i][l+j];
        }
      }
    }
  }
   for(int j=1;j      for(int k=-1;k        for(int l=-1;l          M[0][j]+=m[k+1][l+1]*F[max(k,0)][l+j];
          M[M.dim1()-1][j]+=m[k+1][l+1]*F[M.dim1()-1+min(k,0)][l+j];
        }
      }
  }
  for(int i=1;i      for(int k=-1;k        for(int l=-1;l          M[i][0]+=m[k+1][l+1]*F[k+i][max(l,0)];
          M[i][M.dim2()-1]+=m[k+1][l+1]*F[k+i][M.dim2()-1+min(l,0)];
        }
      }
  }
  for(int k=-1;k    for(int l=-1;l      M[M.dim1()-1][M.dim2()-1]+=m[k+1][l+1]*F[M.dim1()-1+min(k,0)][M.dim2()-1+min(l,0)];
      M[0][M.dim2()-1]+=m[k+1][l+1]*F[max(k,0)][M.dim2()-1+min(l,0)];
      M[M.dim1()-1][0]+=m[k+1][l+1]*F[M.dim1()-1+min(k,0)][max(l,0)];
      M[0][0]+=m[k+1][l+1]*F[max(k,0)][max(l,0)];
    }
  }
  return M;