diff --git a/dataRead.m b/dataRead.m
index 2a4809e..2d3ac79 100644
--- a/dataRead.m
+++ b/dataRead.m
@@ -2,9 +2,10 @@ function [ fsY0, fsY1, fsY2,phaseASpotLoadP,phaseBSpotLoadP,phaseCSpotLoadP ...
     phaseASpotLoadQ,phaseBSpotLoadQ,phaseCSpotLoadQ,setIJ,nodeNum,Balance,phaseABCY] = dataRead(lineZ,dataFile )
 data=dlmread(dataFile);
 baseHighU=data(1,2);
-fsY0=0;
-fsY1=0;
-fsY2=0;
+fsY0=sparse(0);
+fsY1=sparse(0);
+fsY2=sparse(0);
+
 zeroEntries=find(data(:,1)==0);
 lines=data(zeroEntries(1)+1:zeroEntries(2)-1,:);
 [setIJ,nodeNum]=numberNode(lines);
diff --git a/run.m b/run.m
index f50c68d..aa53c83 100644
--- a/run.m
+++ b/run.m
@@ -12,11 +12,6 @@ Tf2p=inv(Tp2f);
 fsY1amp=abs(fsY1);
 [r,c,fsY1ang]=find(fsY1);
 fsY1ang=angle(fsY1ang);
-%%参考电压
-refV0=1;
-refV1=0;
-refV2=0;
-%refV1=1;
 Pabc=phaseASpotLoadP+phaseBSpotLoadP+phaseCSpotLoadP;
 Qabc=phaseASpotLoadQ+phaseBSpotLoadQ+phaseCSpotLoadQ;
 busNum=length(phaseASpotLoadP);
@@ -51,6 +46,9 @@ fsY2=fsY2+sparse(Balance,Balance,ones(length(Balance),1),busNum,busNum);
 fsY0(Balance,:)=0;
 fsY0(:,Balance)=0;
 fsY0=fsY0+sparse(Balance,Balance,ones(length(Balance),1),busNum,busNum);
+%%LU分解
+[fsY0L,fsY0U,fsY0P,fsY0Q,fsY0R]=lu(fsY0);
+[fsY2L,fsY2U,fsY2P,fsY2Q,fsY2R]=lu(fsY2);
 tic
 while(k<=kmax && maxD> EPS)
     k=k+1;
@@ -75,8 +73,11 @@ while(k<=kmax && maxD> EPS)
     If2=conj(f012(3,:)');
     If0(Balance)=0;
     If2(Balance)=0;
-    Vf0=fsY0\If0;
-    Vf2=fsY2\If2;
+    %Vf0=fsY0\If0;
+    %[fsY0L,fsY0U,fsY0P,fsY0Q,fsY0R]
+    Vf0=fsY0Q*(fsY0U\(fsY0L\(fsY0P*(fsY0R\If0))));
+    %Vf2=fsY2\If2;
+    Vf2=fsY2Q*(fsY2U\(fsY2L\(fsY2P*(fsY2R\If2))));
 end
 fprintf('Fortiscue法计算时间 %f\n',toc);
 Vf1=Vmf1.*exp(1j*Vaf1);