基础教程基于吉尔法求解点堆中子动力学方程附MATLAB代码

n=zeros(101,1); C1=zeros(101,1); C2=zeros(101,1); C3=zeros(101,1); C4=zeros(101,1); C5=zeros(101,1); C6=zeros(101,1); n(1)=1; l=0.00002; C1(1)=0.000266/l/0.0127; C2(1)=0.001491/l/0.0317; C3(1)=0.001316/l/0.115; C4(1)=0.002849/l/0.311; C5(1)=0.000896/l/1.4; C6(1)=0.000182/l/3.87; h=0.01; rct=0.003; F=[(rct-0.007)/l,0.0127,0.0317,0.115,0.311,1.4,3.87; 0.000266/l,-0.0127,0,0,0,0,0; 0.001491/l,0,-0.0317,0,0,0,0; 0.001316/l,0,0,-0.115,0,0,0; 0.002849/l,0,0,0,-0.311,0,0; 0.000896/l,0,0,0,0,-1.4,0; 0.000182/l,0,0,0,0,0,-3.87]; A=[1,1,1,1,1,1,1; 0,1,2,3,4,5,6; 0,0,1,3,6,10,15; 0,0,0,1,4,10,20; 0,0,0,0,1,5,15; 0,0,0,0,0,1,6; 0,0,0,0,0,0,1]; I=ones(7); L=[720/1764;1;1624/1764;735/1764;175/1764;21/1764;1/1764]; for i=1:1:100 y=[n(i);C1(i);C2(i);C3(i);C4(i);C5(i);C6(i)]; y1=h.*(F*y); y2=(h^2/factorial(2)).*(F^2*y); y3=(h^3/factorial(3)).*(F^3*y); y4=(h^4/factorial(4)).*(F^4*y); y5=(h^5/factorial(5)).*(F^5*y); y6=(h^6/factorial(6)).*(F^6*y); Z=[y,y1,y2,y3,y4,y5,y6]'; Z0=A*Z; G=-Z0(2,:)'+h.*F*Z0(1,:)'; N=(L(2).*I-(h*L(1)).*F')^(-1); Z1=Z0+dot(L,G').*N; G=-Z1(2,:)'+h.*F*Z1(1,:)'; N=(L(2).*I-(h*L(1)).*F')^(-1); Z2=Z1+dot(L,G').*N; G=-Z2(2,:)'+h.*F*Z2(1,:)'; N=(L(2).*I-(h*L(1)).*F')^(-1); Z3=Z2+dot(L,G').*N; n(i+1)=Z0(1,1); C1(i+1)=Z0(1,2); C2(i+1)=Z0(1,3); C3(i+1)=Z0(1,4); C4(i+1)=Z0(1,5); C5(i+1)=Z0(1,6); C6(i+1)=Z0(1,7); end t=0:0.01:1; plot(t,n)