電力系統(tǒng)穩(wěn)態(tài)實驗報告

1、電力系統(tǒng)穩(wěn)態(tài)潮流計算上機實驗報告一、問題如下圖所示的電力系統(tǒng)網(wǎng)絡(luò)，分別用牛頓拉夫遜法、PQ解耦法、高斯賽德爾法、保留非線性法計算該電力系統(tǒng)的潮流。該電力系統(tǒng)的電路參數(shù)如下，名稱電阻（pu）電抗(pu)-B/2(pu)線路10.040.25-0.25線路20.10.350線路30.080.30-0.25變壓器參數(shù)如下，名稱電阻（pu）電抗(pu)變比變壓器100.031.05:1變壓器200.0151:1.05發(fā)電機的參數(shù)如下，名稱電壓（pu）相角（rad）有功(pu)發(fā)電機11.050*發(fā)電機21.05*5*表示任意值負荷參數(shù)如下，名稱有功(pu)無功(pu)負荷11.60.8負荷22.01.

2、0負荷33.71.3二、問題分析如上圖所示的電力系統(tǒng)，可以看出，節(jié)點1、2、3是PQ節(jié)點，節(jié)點4是PV節(jié)點，而將節(jié)點5作為平衡節(jié)點。根據(jù)問題所需，采用牛頓拉夫遜法、PQ解耦法、高斯賽德爾法、保留非線性法，通過對每次修正量的收斂判據(jù)的判斷，得出整個電力系統(tǒng)的潮流，并分析這四種方法的收斂速度等等。算法分析1.牛頓拉夫遜法節(jié)點5為平衡節(jié)點，不參加整個的迭代過程，節(jié)點1、2、3為PQ節(jié)點，節(jié)點4為PV節(jié)點，計算修正方程中各量，進而得到修正量，判斷修正量是否收斂，如果不收斂，迭代繼續(xù)，如果收斂，算出PQ節(jié)點的電壓幅值以及電壓相角，得出PV節(jié)點的無功量以及電壓相角，得出平衡節(jié)點的輸出功率。潮流方程的直角坐

3、標形式，直角坐標形式的修正方程式，修正方程式中的各量值的計算，Jacobi矩陣的元素計算，牛頓拉夫遜法潮流計算的流程圖如下，2.PQ解耦法如同牛頓拉夫遜法，快速解耦法的前提是，輸電線路的阻抗要比電阻大得多，并且輸電線路兩端的電壓相角相差不大，此時可利用PQ快速解耦法，來計算整個電力系統(tǒng)網(wǎng)絡(luò)的潮流。快速解耦法的迭代方程組，PHQL(U/U)快速解耦法潮流計算的流程圖如下，3.高斯賽德爾法高斯賽德爾法原理較前兩種方法簡單，程序設(shè)計十分容易，占內(nèi)存小，是所有的潮流計算方法中迭代計算量最小的。高斯賽德爾法的迭代格式為，高斯賽德爾法的收斂判據(jù)如下，在高斯賽德爾法中，不應(yīng)對PV節(jié)點的幅值進行修正，只對其電

4、壓相角進行一定的修正。高斯賽德爾法潮流計算的流程圖為，4.保留非線性法保留非線性法主要是在牛頓拉夫遜法的基礎(chǔ)上，通過泰勒展開，保留到二階項，由于三階導(dǎo)數(shù)值等于零，所以泰勒展開式是準確的，無截斷誤差，與牛頓拉夫遜法不同的是，保留非線性法只需算一次jacobi矩陣，每次迭代得到的修正量都是在初始值上的修正量，因此，保留非線性法的計算量小于牛頓拉夫遜法的計算量，大大節(jié)約計算機的內(nèi)存空間，提高計算機的計算速度。保留非線性的迭代格式為，式中，k表示迭代次數(shù)；J為按xx(0)估計而得。收斂判據(jù)為，也可采用相繼二次迭代的二階項之差作為收斂判據(jù)(更合理)，相應(yīng)的收斂判據(jù)如下，保留非線性法的流程圖如下，三、MA

5、TLAB仿真結(jié)果1.牛頓拉夫遜法迭代次數(shù)k=6；各節(jié)點的電壓值、有功功率以及無功功率見下表。名稱電壓幅值電壓相角電壓向量有功功率無功功率節(jié)點10.8683-0.08290.8653 - 0.0719i-1.6-0.8節(jié)點21.07830.31081.0267 + 0.3298i-2-1.0節(jié)點31.0370-0.07461.0341 - 0.0773i-3.7-1.3節(jié)點41.05000.38040.9749 + 0.3899i51.7857節(jié)點51.050001.05002.57602.28032.PQ解耦法迭代次數(shù)k=13；各節(jié)點的電壓值、有功功率以及無功功率見下表。名稱電壓幅值電壓相角電

6、壓向量有功功率無功功率節(jié)點10.8683-0.08290.8653 - 0.0719i-1.6-0.8節(jié)點21.07830.31081.0267 + 0.3298i-2-1.0節(jié)點31.0370-0.07461.0341 - 0.0773i-3.7-1.3節(jié)點41.05000.38040.9749 + 0.3899i51.7857節(jié)點51.050001.05002.57602.28033.高斯賽德爾法迭代次數(shù)k=137;各節(jié)點的電壓值、有功功率以及無功功率見下表。名稱電壓幅值電壓相角電壓向量有功功率無功功率節(jié)點10.8688-0.08470.8657 - 0.0735i-1.6-0.8節(jié)點21

7、.07840.30771.0277 + 0.3266i-2-1.0節(jié)點31.0372-0.07491.0343 - 0.0776i-3.7-1.3節(jié)點41.05000.37720.9762 + 0.3868i51.7857節(jié)點51.050001.05002.57602.28034.保留非線性法 迭代次數(shù)k=11。各節(jié)點的電壓值、有功功率以及無功功率見下表。名稱電壓幅值電壓相角電壓向量有功功率無功功率節(jié)點10.8684-0.08280.8654 - 0.0719i-1.6-0.8節(jié)點21.07830.31081.0267 + 0.3298i-2-1.0節(jié)點31.0370-0.07461.0341

8、 - 0.0773i-3.7-1.3節(jié)點41.05000.38040.9749 + 0.3899i51.7857節(jié)點51.050001.05002.57592.2802四、結(jié)果分析從以上的MATLAB仿真結(jié)果可以看出，牛頓拉夫遜法只需迭代6次，迭代次數(shù)最少，保留非線性迭代11次，大約是牛拉法的兩倍，PQ解耦法迭代次數(shù)13次，收斂速度相比于保留非線性稍慢，而高斯賽德爾法迭代次數(shù)達到137次，高斯賽德爾法算法簡單，占用內(nèi)存小，但是犧牲迭代次數(shù)。從以上的仿真結(jié)果可以得出，牛拉法收斂速度快，算法具有平方收斂特性，是所有算法中收斂最快的，具有良好的收斂可靠性，并且牛頓法所需的內(nèi)存量及每次迭代的時間均較高

9、斯賽德爾法多。在PQ解耦法中，用解兩個階數(shù)幾乎減半的方程組(一個n-1及一個n-m-1)代替牛頓法的結(jié)一個2n-m-2階方程組，顯著地減少了內(nèi)存需求量及計算量，系數(shù)矩陣B及B是兩個常數(shù)陣，為此只需在迭代循環(huán)前一次形成并進行三角分解組成因子表，在迭代過程中反復(fù)應(yīng)用，大大縮短了每次迭代所需時間。快速解耦法達到收斂所需的迭代次數(shù)比牛頓法多，快速解耦法的程序設(shè)計較牛頓法簡單，但從牛頓法到快速解耦法的演化時在元件的R<<X以及線路兩端相角差比較小等假設(shè)基礎(chǔ)上進行的，當系統(tǒng)不符合這些假設(shè)時，迭代就會出現(xiàn)問題。高斯賽德爾法中，原理簡單，程序設(shè)計十分容易，線性非線性方程組均適用，并且導(dǎo)納矩陣是一個

10、對稱且高度稀疏的矩陣，因此占用內(nèi)存非常節(jié)省，每次迭代的計算量也小，是各種潮流算法中最小的。但是收斂速度很慢，迭代次數(shù)將隨所計算網(wǎng)絡(luò)節(jié)點數(shù)的增加而直線上升，從上文的仿真結(jié)果就能看出，收斂速度是四種方法中最慢的。保留非線性法中的雅可比矩陣，只需一次形成，并由三角分解構(gòu)成因子表，而牛頓法中，每次重新形成因子表，保留非線性與牛拉法最大的區(qū)別在于x(k)的含義，在保留非線性中，x(k)是相對于始終不變的初始估計值x(0)的修正量，而在牛拉法中，x(k)是相對于上一次迭代所得到的迭代點x(k)的修正量，但是保留非線性法達到收斂所需迭代次數(shù)多，收斂特性為直線但總計算速度較快。保留非線性法在收斂性方面，屬于“

11、等斜率法”的范疇，和牛頓法的平方收斂特性相比，達到收斂的迭代次數(shù)較牛頓法多，較快速解耦法，收斂的可靠性更好，計算速度可以接近快速解耦法。五、證明 pij=pij+pji證明：因為sij+sji=Ui*Iij*+Uj*Iji* =Ui-Uj*Iij* =Ui-Uj*Ui-UjR+jX* =(Ui-Uj)2R2+X2*(R+jX) =Iij2*R+jX所以pij+pji=Resij+sji=Iij2*R=pij附錄：1. 牛拉法2. clear;3. clc;4. yb=zeros(5,5);5. yb(1,1)=(1.37874-6.26166i)/2;yb(1,2)=-0.62402+3.90

12、015i;yb(1,3)=-0.75471+2.64150i;6. yb(2,2)=(1.45390-66.98082i)/2;yb(2,3)=-0.82987+3.11203i;yb(2,4)=63.49206i;7. yb(3,3)=(1.58459-35.73786i)/2;yb(3,5)=31.74603i;8. yb(4,4)=(-66.66667i)/2;9. yb(5,5)=(-33.33333i)/2;10. yb=yb+conj(yb');11. k=0;12. eps1=10-4;13. jeps=1;14. G=real(yb);15. B=imag(yb);16

13、. e=1; 1 ;1 ;1.05 ;1.05;17. f=zeros(5,1);18. pis=-1.6; -2 ;-3.7; 5;19. qis=-0.8; -1; -1.3;20. deta_p=zeros(4,1);21. deta_q=zeros(3,1);22. deta=zeros(8,1);23. deta_ef=zeros(8,1);24. deta_e=zeros(4,1);25. deta_f=zeros(4,1);26. U=zeros(5,1);27. while (jeps>eps1)28. p=zeros(4,1);29. q=zeros(3,1);30.

14、for i=1:431. for j=1:532. p(i)=p(i)+e(i)*(G(i,j)*e(j)-B(i,j)*f(j)+f(i)*(G(i,j)*f(j)+B(i,j)*e(j);33. end34. deta_p(i)=pis(i)-p(i);35. end36. for i=1:337. for j=1:538. q(i)=q(i)+f(i)*(G(i,j)*e(j)-B(i,j)*f(j)-e(i)*(G(i,j)*f(j)+B(i,j)*e(j);39. end40. deta_q(i)=qis(i)-q(i);41. end42. deta_UU=1.05*1.05-(e

15、(4)*e(4)+f(4)*f(4);43. jacobi=zeros(8,8);44. for i=1:445. for j=1:446. jacobi(2*i-1,2*j-1)=-G(i,j)*e(i)-B(i,j)*f(i);47. jacobi(2*i-1,2*j)=B(i,j)*e(i)-G(i,j)*f(i);48. end49. end50. for i=1:351. for j=1:452. jacobi(2*i,2*j-1)=B(i,j)*e(i)-G(i,j)*f(i);53. jacobi(2*i,2*j)=G(i,j)*e(i)+B(i,j)*f(i);54. end5

16、5. end56. for i=1:2:757. jacobi(i,i)=0;58. jacobi(i,i+1)=0;59. for j=1:560. jacobi(i,i)=jacobi(i,i)-(G(i+1)/2,j)*e(j)-B(i+1)/2,j)*f(j);61. jacobi(i,i+1)=jacobi(i,i+1)-(G(i+1)/2,j)*f(j)+B(i+1)/2,j)*e(j);62. end63. jacobi(i,i)=jacobi(i,i)-G(i+1)/2,(i+1)/2)*e(i+1)/2)-B(i+1)/2,(i+1)/2)*f(i+1)/2);64. jac

17、obi(i,i+1)=jacobi(i,i+1)+B(i+1)/2,(i+1)/2)*e(i+1)/2)-G(i+1)/2,(i+1)/2)*f(i+1)/2);65. end66. for i=2:2:667. jacobi(i,i-1)=0;68. jacobi(i,i)=0;69. for j=1:570. jacobi(i,i-1)=jacobi(i,i-1)-(G(i/2,j)*f(j)-B(i/2,j)*e(j);71. jacobi(i,i)=jacobi(i,i)-(G(i/2,j)*e(j)-B(i/2,j)*f(j);72. end73. jacobi(i,i-1)=jac

18、obi(i,i-1)+B(i/2,i/2)*e(i/2)-G(i/2,i/2)*f(i/2);74. jacobi(i,i)=jacobi(i,i)+G(i/2,i/2)*e(i/2)+B(i/2,i/2)*f(i/2);75. end76. jacobi(8,7)=-2*e(4);77. jacobi(8,8)=-2*f(4);78. for i=1:2:779. deta(i)=deta_p(i+1)/2);80. end81. for i=2:2:682. deta(i)=deta_q(i/2);83. end84. deta(8)=deta_UU;85. deta_ef=-inv(ja

19、cobi)*deta;86. for i=1:487. deta_e(i)=deta_ef(2*i-1);88. deta_f(i)=deta_ef(2*i);89. end90. for i=1:491. e(i)=e(i)+deta_e(i);92. f(i)=f(i)+deta_f(i);93. end94. jeps=max(max(abs(deta_ef),max(abs(deta_p);95. jeps=max(jeps,max(abs(deta_q);96. k=k+1;97. end98. for i=1:599. U(i)=e(i)+sqrt(-1)*f(i);100. en

20、d101. disp('µçѹµÄ·ùÖµ£º');102. disp(abs(U);103. disp('µçѹ£º');104. disp(U);105. theta=zeros(5,1);106. for i=1:5107. theta(i)=angle(U(i);108. end109. disp(theta);110. y0=-4i;-8i;-4i;0;0;111. s=

21、zeros(5,5);112. for i=1:5113. for j=1:5114. s(i,j)=U(i)*(conj(U(i)*conj(y0(i)+(conj(U(i)-conj(U(j)*conj(-yb(i,j);115. end116. end117. for i=1:5118. s(i,i)=0;119. end120. deta_s=zeros(5,5);121. p=zeros(5,1);122. q=zeros(5,1);123. for i=1:5124. for j=1:5125. p(i)=p(i)+e(i)*(G(i,j)*e(j)-B(i,j)*f(j)+f(i

22、)*(G(i,j)*f(j)+B(i,j)*e(j);126. end127. end128. for i=1:5129. for j=1:5130. q(i)=q(i)+f(i)*(G(i,j)*e(j)-B(i,j)*f(j)-e(i)*(G(i,j)*f(j)+B(i,j)*e(j);131. end132. end133. for i=1:5134. for j=1:5135. deta_s(i,j)=s(i,j)+s(j,i);136. end 137. end138. disp(s);139. disp(deta_s);140. disp(p);disp(q); 2.PQ解耦法cl

23、ear;clc;yb=zeros(5,5);yb(1,1)=(1.37874-6.26166i)/2;yb(1,2)=-0.62402+3.90015i;yb(1,3)=-0.75471+2.64150i;yb(2,2)=(1.45390-66.98082i)/2;yb(2,3)=-0.82987+3.11203i;yb(2,4)=63.49206i;yb(3,3)=(1.58459-35.73786i)/2;yb(3,5)=31.74603i;yb(4,4)=(-66.66667i)/2;yb(5,5)=(-33.33333i)/2;yb=yb+conj(yb');k=1;kp=0;

24、kq=0;eps1=10-5;G=real(yb);B=imag(yb);B1=B(1:4,1:4);B2=B(1:3,1:3);U=1; 1 ;1 ;1.05 ;1.05;theta=zeros(5,1);deta_theta=zeros(4,1);deta_U=zeros(3,1);pis=-1.6; -2 ;-3.7; 5;qis=-0.8; -1; -1.3;deta_p=zeros(4,1);deta_q=zeros(3,1);while (kp=0|kq=0) p=zeros(4,1); for i=1:4 for j=1:5 p(i)=p(i)+U(i)*U(j)*(G(i,j)

25、*cos(theta(i)-theta(j)+B(i,j)*sin(theta(i)-theta(j); end deta_p(i)=pis(i)-p(i); end for i=1:4 deta_p(i)=deta_p(i)/U(i); end deta_theta=-inv(B1)*deta_p; for i=1:4 deta_theta(i)=deta_theta(i)/U(i); theta(i)=theta(i)+deta_theta(i); end eps_theta=max(max(abs(deta_theta),max(abs(deta_p); if eps_theta<

26、eps1 kp=1; else kp=0; end q=zeros(3,1); for i=1:3 for j=1:5 q(i)=q(i)+U(i)*U(j)*(G(i,j)*sin(theta(i)-theta(j)-B(i,j)*cos(theta(i)-theta(j); end deta_q(i)=qis(i)-q(i); end for i=1:3 deta_q(i)=deta_q(i)/U(i); end deta_U=-inv(B2)*deta_q; for i=1:3 U(i)=U(i)+deta_U(i); end eps_U=max(max(abs(deta_U),max(

27、abs(deta_q); if eps_U<eps1 kq=1; else kq=1; end if max(eps_theta,eps_U)>eps1 k=k+1; endendfor i=1:5 U(i)=U(i)*(cos(theta(i)+sqrt(-1)*sin(theta(i);enddisp(U);disp(theta);p=zeros(5,1);q=zeros(5,1);for i=1:5 for j=1:5 p(i)=p(i)+abs(U(i)*abs(U(j)*(G(i,j)*cos(theta(i)-theta(j)+B(i,j)*sin(theta(i)-t

28、heta(j); endenddisp(p);for i=1:5 for j=1:5 q(i)=q(i)+abs(U(i)*abs(U(j)*(G(i,j)*sin(theta(i)-theta(j)-B(i,j)*cos(theta(i)-theta(j); endenddisp(q);y0=-4i;-8i;-4i;0;0;s=zeros(5,5);for i=1:5 for j=1:5 s(i,j)=U(i)*(conj(U(i)*conj(y0(i)+(conj(U(i)-conj(U(j)*conj(-yb(i,j); endendfor i=1:5 s(i,i)=0;enddeta_

29、s=zeros(5,5);for i=1:5 for j=1:5 deta_s(i,j)=s(i,j)+s(j,i); end enddisp(s);disp(deta_s);3.高斯賽德爾法clear;clc;yb=zeros(5,5);yb(1,1)=(1.37874-6.26166i)/2;yb(1,2)=-0.62402+3.90015i;yb(1,3)=-0.75471+2.64150i;yb(2,2)=(1.45390-66.98082i)/2;yb(2,3)=-0.82987+3.11203i;yb(2,4)=63.49206i;yb(3,3)=(1.58459-35.73786

30、i)/2;yb(3,5)=31.74603i;yb(4,4)=(-66.66667i)/2;yb(5,5)=(-33.33333i)/2;yb=yb+conj(yb');G=real(yb);B=imag(yb);s=-1.6-0.8i;-2-i;-3.7-1.3i;5;U0=1;1;1;1.05;1.05;U1=zeros(4,1);deta_U=zeros(4,1);eps1=10-4;eps_U=1;theta=0;k=0;while (eps_U>eps1) q=0; for j=1:5 q=q+U0(4)*conj(yb(4,j)*conj(U0(j); end q=i

31、mag(q); s(4)=5+sqrt(-1)*q; p=0; for j=1:5 p=p+yb(4,j)*U0(j); end p=p-yb(4,4)*U0(4); U1(4)=(conj(s(4)/conj(U0(4)-p)/yb(4,4); theta=angle(U1(4); U1(4)=1.05*(cos(theta)+sqrt(-1)*sin(theta); for i=1:3 a=0; b=0; for j=1:i-1 a=a+yb(i,j)*U1(j); end for j=i+1:4 b=b+yb(i,j)*U0(j); end U1(i)=(conj(s(i)/conj(U

32、0(i)-yb(i,5)*U0(5)-a-b)/yb(i,i); end for i=1:4 deta_U=U1(i)-U0(i); U0(i)=U1(i); end eps_U=max(abs(deta_U); k=k+1;enddisp(abs(U0);disp(U0);theta=zeros(5,1);for i=1:5 theta(i)=angle(U0(i);enddisp(theta);y0=-4i;-8i;-4i;0;0;s=zeros(5,5);for i=1:5 for j=1:5 s(i,j)=U0(i)*(conj(U0(i)*conj(y0(i)+(conj(U0(i)

33、-conj(U0(j)*conj(-yb(i,j); endendfor i=1:5 s(i,i)=0;enddeta_s=zeros(5,5);for i=1:5 for j=1:5 deta_s(i,j)=s(i,j)+s(j,i); end enddisp(s);disp(deta_s);4.保留非線性法clear;clc;yb=zeros(5,5);yb(1,1)=(1.37874-6.26166i)/2;yb(1,2)=-0.62402+3.90015i;yb(1,3)=-0.75471+2.64150i;yb(2,2)=(1.45390-66.98082i)/2;yb(2,3)=-

34、0.82987+3.11203i;yb(2,4)=63.49206i;yb(3,3)=(1.58459-35.73786i)/2;yb(3,5)=31.74603i;yb(4,4)=(-66.66667i)/2;yb(5,5)=(-33.33333i)/2;yb=yb+conj(yb');k=0;eps1=10-4;unlineeps=1;G=real(yb);B=imag(yb);U=zeros(5,1);e=1; 1 ;1 ;1.05 ;1.05;f=zeros(5,1);pis=-1.6; -2 ;-3.7; 5;qis=-0.8; -1; -1.3;deta_p=zeros(4

35、,1);deta_q=zeros(3,1);deta=zeros(8,1);deta_ef=zeros(8,1);deta_e0=zeros(5,1);deta_e1=zeros(5,1);deta_f0=zeros(5,1);deta_f1=zeros(5,1);p=zeros(4,1);q=zeros(3,1);for i=1:4 for j=1:5 p(i)=p(i)+e(i)*(G(i,j)*e(j)-B(i,j)*f(j)+f(i)*(G(i,j)*f(j)+B(i,j)*e(j); end deta_p(i)=pis(i)-p(i);endfor i=1:3 for j=1:5 q

36、(i)=q(i)+f(i)*(G(i,j)*e(j)-B(i,j)*f(j)-e(i)*(G(i,j)*f(j)+B(i,j)*e(j); end deta_q(i)=qis(i)-q(i);enddeta_UU=1.05*1.05-(e(4)*e(4)+f(4)*f(4);jacobi=zeros(8,8);for i=1:4 for j=1:4 jacobi(2*i-1,2*j-1)=-G(i,j)*e(i)-B(i,j)*f(i); jacobi(2*i-1,2*j)=B(i,j)*e(i)-G(i,j)*f(i); endendfor i=1:3 for j=1:4 jacobi(2*

37、i,2*j-1)=B(i,j)*e(i)-G(i,j)*f(i); jacobi(2*i,2*j)=G(i,j)*e(i)+B(i,j)*f(i); endendfor i=1:2:7 jacobi(i,i)=0; jacobi(i,i+1)=0; for j=1:5 jacobi(i,i)=jacobi(i,i)-(G(i+1)/2,j)*e(j)-B(i+1)/2,j)*f(j); jacobi(i,i+1)=jacobi(i,i+1)-(G(i+1)/2,j)*f(j)+B(i+1)/2,j)*e(j); end jacobi(i,i)=jacobi(i,i)-G(i+1)/2,(i+1

38、)/2)*e(i+1)/2)-B(i+1)/2,(i+1)/2)*f(i+1)/2); jacobi(i,i+1)=jacobi(i,i+1)+B(i+1)/2,(i+1)/2)*e(i+1)/2)-G(i+1)/2,(i+1)/2)*f(i+1)/2);endfor i=2:2:6 jacobi(i,i-1)=0; jacobi(i,i)=0; for j=1:5 jacobi(i,i-1)=jacobi(i,i-1)-(G(i/2,j)*f(j)-B(i/2,j)*e(j); jacobi(i,i)=jacobi(i,i)-(G(i/2,j)*e(j)-B(i/2,j)*f(j); end

39、 jacobi(i,i-1)=jacobi(i,i-1)+B(i/2,i/2)*e(i/2)-G(i/2,i/2)*f(i/2); jacobi(i,i)=jacobi(i,i)+G(i/2,i/2)*e(i/2)+B(i/2,i/2)*f(i/2);endjacobi(8,7)=-2*e(4);jacobi(8,8)=-2*f(4);while (unlineeps>eps1) deta_pp=zeros(4,1); for i=1:4 for j=1:5 deta_pp(i)=deta_pp(i)+deta_e0(i)*(G(i,j)*deta_e0(j)-B(i,j)*deta_f

40、0(j)+deta_f0(i)*(G(i,j)*deta_f0(j)+B(i,j)*deta_e0(j); end end deta_qq=zeros(3,1); for i=1:3 for j=1:5 deta_qq(i)=deta_qq(i)+deta_f0(i)*(G(i,j)*deta_e0(j)-B(i,j)*deta_f0(j)-deta_e0(i)*(G(i,j)*deta_f0(j)+B(i,j)*deta_e0(j); end end deta_UUU=-(deta_e0(4)*deta_e0(4)+deta_f0(4)*deta_f0(4); for i=1:2:7 deta(i)=deta_p(i+1)/2)-deta_pp(i+1)/2); end for i=2:2:6 deta(i)=deta_q(i/2)-deta_qq(i/2); end deta(