LECCION 9 - SOLUCION ECUACIONES

Off[General::"spell1"]
Off[General::"spell"]

01 - FACTORIZACION DE UNA SKYMATRIZ SIMETRICA

SymmSkyMatrixFactor[S_,tol_]:=Module[{p,a,fail,i,j,k,l,m,n,ii,ij,jj,jk,jmj,d,s,row,v},row=SymmSkyMatrixRowLengths[S];s=Max[row];{p,a}=S;n=Length[p]-1;v=Table[0,{n}];fail=0;Do[jj=p[[j+1]];If[jj<0|row[[j]]==0,Continue[]];d=a[[jj]];jmj=Abs[p[[j]]];jk=jj-jmj;Do[i=j-jk+k;v[[k]]=0;ii=p[[i+1]];If[ii<0,Continue[]];m=Min[ii-Abs[p[[i]]],k]-1;ij=jmj+k;v[[k]]=a[[ij]];v[[k]]-=Take[a,{ii-m,ii-1}].Take[v,{k-m,k-1}];a[[ij]]=v[[k]]*a[[ii]],{k,1,jk-1}];d-=Take[a,{jmj+1,jmj+jk-1}].Take[v,{1,jk-1}];If[Abs[d]<tol*row[[j]],fail=j;a[[jj]]=Infinity;Break[]];a[[jj]]=1/d,{j,1,n}];Return[{{p,a},fail}]];

SymmSkyMatrixRowLengths[S_]:=Module[{p,a,i,j,n,ii,jj,m,d,row},{p,a}=S;n=Length[p]-1;row=Table[0,{n}];Do[ii=p[[i+1]];If[ii<0,Continue[]];m=ii-i;row[[i]]=a[[ii]]^2;Do[If[p[[j+1]]>0,d=a[[m+j]]^2;row[[i]]+=d;row[[j]]+=d],{j,Max[1,Abs[p[[i]]]-m+1],Min[n,i]-1}],{i,1,n}];Return[Sqrt[row]];];

02a - SOLUCION DEL SISTEMA DE ECUACIONES A PARTIR DE LA SKYMATRIZ FACTORIZADA Y DEL VECTOR DE TERMINOS INDEPENDIENTES

SymmSkyMatrixVectorSolve[S_,b_]:=Module[{p,a,n,i,j,k,m,ii,jj,bi,x},{p,a}=S;n=Length[p]-1;x=b;If[n≠Length[x],Print["Inconsistentmatrixdimensionsin","SymmSkyMatrixVectorSolve"];Return[Null]];Do[ii=p[[i+1]];If[ii≥0,Continue[]];ii=-ii;k=i-ii+Abs[p[[i]]]+1;bi=x[[i]];If[bi==0,Continue[]];Do[jj=p[[j+1]],If[jj<0,Continue[]];m=j-i;If[m<0,x[[j]]-=a[[ii+m]]*bi;Break[]];ij=jj-m;If[ij>Abs[p[[j]]],x[[j]]-=a[[ij]]*bi],{j,k,n}],{i,1,n}];Do[ii=p[[i+1]];If[ii<0,x[[i]]=0;Continue[]];imi=Abs[p[[i]]];m=ii-imi-1;x[[i]]-=Take[a,{imi+1,imi+m}].Take[x,{i-m,i-1}],{i,1,n}];Do[ii=Abs[p[[i+1]]];x[[i]]*=a[[ii]],{i,1,n}];Do[ii=p[[i+1]];If[ii<0,x[[i]]=b[[i]];Continue[]];m=ii-Abs[p[[i]]]-1;Do[x[[i-j]]-=a[[ii-j]]*x[[i]],{j,1,m}],{i,n,1,-1}];Return[x]];

02b -SOLUCION PARA TODOS LOS VECTORES DE CARGA CORRESPONDIENTES A LOS DIFERENTES CASOS DE CARGA

SymmSkyMatrixColBlockSolve[S_,b_]:=Module[{p,a,n,nrhs,i,j,k,m,r,ii,jj,bi,x},{p,a}=S;n=Length[p]-1;x=b;If[n≠Dimensions[x][[1]],Print["Inconsistentmatrixdimensionsin","SymmSkyMatrixBlockColSolve"];Return[Null]];nrhs=Dimensions[x][[2]];Do[ii=p[[i+1]];If[ii≥0,Continue[]];ii=-ii;k=i-ii+Abs[p[[i]]]+1;Do[bi=x[[i,r]];If[bi==0,Continue[]];Do[jj=p[[j+1]],If[jj<0,Continue[]];m=j-i;If[m<0,x[[j,r]]-=a[[ii+m]]*bi;Break[]];ij=jj-m;If[ij>Abs[p[[j]]],x[[j,r]]-=a[[ij]]*bi],{j,k,n}],{r,1,nrhs}],{i,1,n}];Do[ii=p[[i+1]];If[ii<0,Do[x[[i,r]]=0,{r,1,nrhs}];Continue[]];imi=Abs[p[[i]]];m=ii-imi-1;Do[Do[x[[i,r]]-=a[[imi+j]]*x[[i-m+j-1,r]],{j,1,m}],{r,1,nrhs}],{i,1,n}];Do[ii=Abs[p[[i+1]]];Do[x[[i,r]]*=a[[ii]],{r,1,nrhs}],{i,1,n}];Do[ii=p[[i+1]];If[ii<0,Do[x[[i,r]]=b[[i,r]],{r,1,nrhs}];Continue[]];m=ii-Abs[p[[i]]]-1;Do[Do[x[[i-j,r]]-=a[[ii-j]]*x[[i,r]],{j,1,m}],{r,1,nrhs}],{i,n,1,-1}];Return[x]];

SymmSkyMatrixRowBlockSolve[S_,b_]:=Module[{p,a,n,nrhs,i,j,k,m,r,ii,jj,bi,x},{p,a}=S;n=Length[p]-1;x=b;If[n≠Dimensions[x][[2]],Print["Inconsistentmatrixdimensionsin","SymmSkyMatrixBlockRowSolve"];Return[Null]];nrhs=Dimensions[x][[1]];Do[ii=p[[i+1]];If[ii≥0,Continue[]];ii=-ii;k=i-ii+Abs[p[[i]]]+1;Do[bi=x[[r,i]];If[bi==0,Continue[]];Do[jj=p[[j+1]],If[jj<0,Continue[]];m=j-i;If[m<0,x[[j,r]]-=a[[ii+m]]*bi;Break[]];ij=jj-m;If[ij>Abs[p[[j]]],x[[r,j]]-=a[[ij]]*bi],{j,k,n}],{r,1,nrhs}],{i,1,n}];Do[ii=p[[i+1]];If[ii<0,Do[x[[r,i]]=0,{r,1,nrhs}];Continue[]];imi=Abs[p[[i]]];m=ii-imi-1;Do[Do[x[[r,i]]-=a[[imi+j]]*x[[r,i-m+j-1]],{j,1,m}],{r,1,nrhs}],{i,1,n}];Do[ii=Abs[p[[i+1]]];Do[x[[r,i]]*=a[[ii]],{r,1,nrhs}],{i,1,n}];Do[ii=p[[i+1]];If[ii<0,Do[x[[r,i]]=b[[r,i]],{r,1,nrhs}];Continue[]];m=ii-Abs[p[[i]]]-1;Do[Do[x[[r,i-j]]-=a[[ii-j]]*x[[r,i]],{j,1,m}],{r,1,nrhs}],{i,n,1,-1}];Return[x]];

03a -  MULTIPLICACION DE UNA SKYMATRIZ POR VECTOR DESPLAZAMIENTOS NODALES PARA OBTENCION FUERZAS NODALES RESULTANTES

SymmSkyMatrixVectorMultiply[S_,x_]:=Module[{p,a,n,i,j,k,m,ii,b},{p,a}=S;n=Length[p]-1;If[n≠Length[x],Print["Inconsistentmatrixdimensionsin","SymmSkyMatrixVectorMultiply"];Return[Null]];b=Table[a[[Abs[p[[i+1]]]]]*x[[i]],{i,1,n}];Do[ii=Abs[p[[i+1]]];m=ii-Abs[p[[i]]]-1;If[m≤0,Continue[]];b[[i]]+=Take[a,{ii-m,ii-1}].Take[x,{i-m,i-1}];Do[b[[i-k]]+=a[[ii-k]]*x[[i]],{k,1,m}],{i,1,n}];Return[b]];

03b -  MULTIPLICACION DE UNA SKYMATRIZ POR SOLUCIONES VECTORES DESPLAZAMIENTOS PARA VARIOS CASOS DE CARGA - OBTENCION VECTORES FUERZAS NODALES RESULTANTES.

SymmSkyMatrixColBlockMultiply[S_,x_]:=Module[{p,a,n,nrhs,i,j,k,m,r,ii,aij,b},{p,a}=S;n=Length[p]-1;If[n≠Dimensions[x][[1]],Print["Inconsistentmatrixdimensionsin","SymmSkyMatrixColBlockMultiply"];Return[Null]];nrhs=Dimensions[x][[2]];b=Table[0,{n},{nrhs}];Do[ii=Abs[p[[i+1]]];m=ii-Abs[p[[i]]]-1;Do[b[[i,r]]=a[[ii]]*x[[i,r]],{r,1,nrhs}];Do[j=i-k;aij=a[[ii-k]];If[aij==0,Continue[]];Do[b[[i,r]]+=aij*x[[j,r]];b[[j,r]]+=aij*x[[i,r]],{r,1,nrhs}],{k,1,m}],{i,1,n}];Return[b]];

SymmSkyMatrixRowBlockMultiply[S_,x_]:=Module[{p,a,n,nrhs,i,j,k,m,r,ii,aij,b},{p,a}=S;n=Length[p]-1;If[n≠Dimensions[x][[2]],Print["Inconsistentmatrixdimensionsin","SymmSkyMatrixRowBlockMultiply"];Return[Null]];nrhs=Dimensions[x][[1]];b=Table[0,{nrhs},{n}];Do[ii=Abs[p[[i+1]]];m=ii-Abs[p[[i]]]-1;Do[b[[r,i]]=a[[ii]]*x[[r,i]],{r,1,nrhs}];Do[j=i-k;aij=a[[ii-k]];If[aij==0,Continue[]];Do[b[[r,i]]+=aij*x[[r,j]];b[[r,j]]+=aij*x[[r,i]],{r,1,nrhs}],{k,1,m}],{i,1,n}];Return[b]];

04 - VISUALIZACION DE SKYMATRICES

SymmSkyMatrixLowerTrianglePrint[S_]:=Module[{p,a,cycle,i,ii,ij,it,j,jj,j1,j2,jref,jbeg,jend,jt,kcmax,kc,kr,m,n,c,t},{p,a}=S;n=Dimensions[p][[1]]-1;kcmax=5;jref=0;Label[cycle];Print[""];jbeg=jref+1;jend=Min[jref+kcmax,n];kc=jend-jref;t=Table["",{n-jref+1},{kc+1}];Do[If[p[[j+1]]>0,c="",c="*"];t[[1,j-jref+1]]=StringJoin[c,"Col",ToString[PaddedForm[j,3]]],{j,jbeg,jend}];it=1;Do[ii=Abs[p[[i+1]]];m=ii-Abs[p[[i]]]-1;j1=Max[i-m,jbeg];j2=Min[i,jend];kr=j2-j1+1;If[kr≤0,Continue[]];If[p[[i+1]]>0,c="",c="*"];it++;t[[it,1]]=StringJoin[c,"Row",ToString[PaddedForm[i,3]]];jt=j1-jbeg+2;ij=j1+ii-i;Do[t[[it,jt++]]=PaddedForm[a[[ij++]]//FortranForm,{7,4}],{j,1,kr}],{i,jbeg,n}];Print[TableForm[Take[t,it],TableAlignments->{Right,Right},TableDirections->{Column,Row},TableSpacing->{0,2}]];jref=jend;If[jref<n,Goto[cycle]];];

SymmSkyMatrixUpperTrianglePrint[S_]:=Module[{p,a,cycle,i,ij,it,j,j1,j2,jref,jbeg,jend,kcmax,k,kc,m,n,c,t},{p,a}=S;n=Dimensions[p][[1]]-1;kcmax=5;jref=0;Label[cycle];Print[""];jbeg=jref+1;jend=Min[jref+kcmax,n];kc=jend-jref;t=Table["",{jend+1},{kc+1}];Do[If[p[[j+1]]>0,c="",c="*"];t[[1,j-jref+1]]=StringJoin[c,"Col",ToString[PaddedForm[j,3]]],{j,jbeg,jend}];it=1;Do[it++;If[p[[i+1]]>0,c="",c="*"];t[[it,1]]=StringJoin[c,"Row",ToString[PaddedForm[i,3]]];j=jref;Do[j++;If[j<i,Continue[]];ij=Abs[p[[j+1]]]+i-j;If[ij≤Abs[p[[j]]],Continue[]];t[[it,k+1]]=PaddedForm[a[[ij]]//FortranForm,{7,4}],{k,1,kc}],{i,1,jend}];Print[TableForm[Take[t,it],TableAlignments->{Right,Right},TableDirections->{Column,Row},TableSpacing->{0,2}]];jref=jend;If[jref<n,Goto[cycle]];];

05 -  VISUALIZACION TERMINOS NO NULOS DE SKYMATRICES - TERMINOS A ALMACENAR EN VECTOR S

SymmSkyMatrixLowerTriangleMap[S_]:=Module[{p,a,cycle,i,ii,ij,it,itop,j,jj,j1,j2,jref,jbeg,jend,jt,kcmax,kc,kr,m,n,c,t},{p,a}=S;n=Dimensions[p][[1]]-1;kcmax=40;jref=0;Label[cycle];Print[""];jbeg=jref+1;jend=Min[jref+kcmax,n];kc=jend-jref;itop=2;If[jend>9,itop=3];If[jend>99,itop=4];If[jend>999,itop=5];t=Table["",{n-jref+itop},{kc+1}];it=0;If[itop≥5,it++;Do[m=Floor[j/1000];If[m>0,t[[it,j-jref+1]]=ToString[Mod[m,10]]],{j,jbeg,jend}]];If[itop≥4,it++;Do[m=Floor[j/100];If[m>0,t[[it,j-jref+1]]=ToString[Mod[m,10]]],{j,jbeg,jend}]];If[itop≥3,it++;Do[m=Floor[j/10];If[m>0,t[[it,j-jref+1]]=ToString[Mod[m,10]]],{j,jbeg,jend}]];it++;Do[t[[it,j-jref+1]]=ToString[Mod[j,10]],{j,jbeg,jend}];it++;Do[If[p[[j+1]]<0,t[[it,j-jref+1]]="*"],{j,jbeg,jend}];Do[ii=Abs[p[[i+1]]];m=ii-Abs[p[[i]]]-1;j1=Max[i-m,jbeg];j2=Min[i,jend];kr=j2-j1+1;If[kr≤0,Continue[]];If[p[[i+1]]>0,c="",c="*"];it++;t[[it,1]]=StringJoin[ToString[PaddedForm[i,2]],c];jt=j1-jbeg+2;ij=j1+ii-i;Do[c="0";If[a[[ij]]>0,c="+"];If[a[[ij++]]<0,c="-"];t[[it,jt++]]=c,{j,1,kr}],{i,jbeg,n}];Print[TableForm[Take[t,it],TableAlignments->{Right,Right},TableDirections->{Column,Row},TableSpacing->{0,0}]];jref=jend;If[jref<n,Goto[cycle]];];

SymmSkyMatrixUpperTriangleMap[S_]:=Module[{p,a,cycle,i,ij,it,itop,j,j1,j2,jref,jbeg,jend,kcmax,k,kc,m,n,c,t},{p,a}=S;n=Dimensions[p][[1]]-1;kcmax=40;jref=0;Label[cycle];Print[""];jbeg=jref+1;jend=Min[jref+kcmax,n];kc=jend-jref;itop=2;If[jend>9,itop=3];If[jend>99,itop=4];If[jend>999,itop=5];t=Table["",{jend+itop},{kc+1}];it=0;If[itop≥5,it++;Do[m=Floor[j/1000];If[m>0,t[[it,j-jref+1]]=ToString[Mod[m,10]]],{j,jbeg,jend}]];If[itop≥4,it++;Do[m=Floor[j/100];If[m>0,t[[it,j-jref+1]]=ToString[Mod[m,10]]],{j,jbeg,jend}]];If[itop≥3,it++;Do[m=Floor[j/10];If[m>0,t[[it,j-jref+1]]=ToString[Mod[m,10]]],{j,jbeg,jend}]];it++;Do[t[[it,j-jref+1]]=ToString[Mod[j,10]],{j,jbeg,jend}];it++;Do[If[p[[j+1]]<0,t[[it,j-jref+1]]="*"],{j,jbeg,jend}];Do[it++;If[p[[i+1]]>0,c="",c="*"];t[[it,1]]=StringJoin[ToString[PaddedForm[i,2]],c];j=jref;Do[j++;If[j<i,Continue[]];ij=Abs[p[[j+1]]]+i-j;If[ij≤Abs[p[[j]]],Continue[]];c="0";If[a[[ij]]>0,c="+"];If[a[[ij++]]<0,c="-"];t[[it,k+1]]=c,{k,1,kc}],{i,1,jend}];Print[TableForm[Take[t,it],TableAlignments->{Right,Right},TableDirections->{Column,Row},TableSpacing->{0,0}]];jref=jend;If[jref<n,Goto[cycle]];];

06 -  RECONSTRUCCION DE SKYMATRICES A PARTIR DE FACTORES

SymmSkyMatrixLDUReconstruct[S_]:=Module[{p,ldu,a,v,n,i,ii,ij,j,jj,jk,jmj,k,m},{p,ldu}=S;a=ldu;n=Length[p]-1;v=Table[0,{n}];Do[jmj=Abs[p[[j]]];jj=p[[j+1]];If[jj<0,Continue[]];jk=jj-jmj;v[[jk]]=ldu[[jj]];Do[ij=jmj+k;i=j+ij-jj;ii=p[[i+1]];If[ii<0,v[[k]]=0;Continue[]];If[i≠j,v[[k]]=ldu[[ij]]*ldu[[ii]]];m=Min[ii-Abs[p[[i]]],k];a[[ij]]=v[[k]];a[[ij]]+=Take[ldu,{ii-m+1,ii-1}].Take[v,{k-m+1,k-1}],{k,1,jk}],{j,1,n}];Return[{p,a}];];

SymmSkyMatrixLDinvUReconstruct[S_]:=Module[{p,ldu,a,v,n,i,ii,ij,j,jj,jk,jmj,k,m},{p,ldu}=S;a=ldu;n=Length[p]-1;v=Table[0,{n}];Do[jmj=Abs[p[[j]]];jj=p[[j+1]];If[jj<0,Continue[]];jk=jj-jmj;v[[jk]]=1/ldu[[jj]];Do[ij=jmj+k;i=j+ij-jj;ii=p[[i+1]];If[ii<0,v[[k]]=0;Continue[]];If[i≠j,v[[k]]=ldu[[ij]]/ldu[[ii]]];m=Min[ii-Abs[p[[i]]],k];a[[ij]]=v[[k]];a[[ij]]+=Take[ldu,{ii-m+1,ii-1}].Take[v,{k-m+1,k-1}],{k,1,jk}],{j,1,n}];Return[{p,a}];];

07 -  UTILIDADES ADICIONALES

SymmSkyMatrixConvertToFull[S_]:=Module[{p,a,aa,n,j,jj,jmj,k},{p,a}=S;n=Length[p]-1;aa=Table[0,{n},{n}];Do[jmj=Abs[p[[j]]];jj=Abs[p[[j+1]]];aa[[j,j]]=a[[jj]];Do[aa[[j,j-k]]=aa[[j-k,j]]=a[[jj-k]],{k,1,jj-jmj-1}],{j,1,n}];Return[aa];];

SymmSkyMatrixConvertUnitUpperTriangleToFull[S_]:=Module[{p,ldu,aa,n,j,jj,jmj,k},{p,ldu}=S;n=Length[p]-1;aa=Table[0,{n},{n}];Do[jmj=Abs[p[[j]]];jj=Abs[p[[j+1]]];aa[[j,j]]=1;Do[aa[[j-k,j]]=ldu[[jj-k]],{k,1,jj-jmj-1}],{j,1,n}];Return[aa];];

SymmSkyMatrixConvertDiagonalToFull[S_]:=Module[{p,ldu,aa,n,i,j,jj,jmj,k},{p,ldu}=S;n=Length[p]-1;aa=Table[0,{n},{n}];Do[jj=Abs[p[[j+1]]];aa[[j,j]]=ldu[[jj]],{j,1,n}];Return[aa];];