正交配置求解问题：运用正交配置法求解有轴向扩散的固定床反应器中催化反应的温度和浓度分布。柱形固体床反应器中催化反应的温度和浓度方程为：＝+=+︱r=1=|r=0=0-|r=1=Biw[T(1,z)-Tw(z)],-|r=1=0T(r,0)=T0,c(r,0)=c0＝+=+︱r=1=|r=0=0-|r=1=Biw[T(1,z)-Tw(z)],-|r=1=0T(r,0)=T0,c(r,0)=c0其中R（c，T）为催化反应的速率方程，其形式为R(c,T)=+=+︱r=1=|r=0=0-|r=1=Biw[T(1,z)-Tw(z)],-|r=1=0T(r,0)=T0,c(r,0)=c0其中R（c，T）为催化反应的速率方程，其形式为R(c,T)=解题思路：应用对称的正交配置法，有下面的方程和初始条件：=+(1-)=+(1-)Tj(0)=T0,cj(0)=c0边界条件为：-AN+1,iTi=Biw(TN+1-Tw),AN+1,ici=0将温度和浓度的边界条件代入微分方程，消去边界值，可得2N个常微分方程，而将两边界条件的代数方程同2N个常微分方程组联合，就组成2N+2个微分代数方程组。结合正交配置系数的计算程序与常微分方程组或微分方程组求解程序，可得到反应器中的温度和浓度分布。具体做法如下：一、利用对称的正交配置格式：1、对称常微分方程程序：（COLLAB.FOR,DLSODE.FOR）主程序：IMPLICITREAL*8(A-H,O-Z)EXTERNALFEX,JEXDIMENSIONAS(19,19),BS(19,19),Q(19,19),XS(19),WS(19)DIMENSIONDIF1(19),DIF2(19),DIF3(19),ROOT(19),V1(19),V2(19)DIMENSIONY(99),ATOL(99),RWORK(10920),IWORK(120)DOUBLEPRECISIONYN1,YN2COMMON/AB/N,AS,BSCOMMON/BC/YN1,YN2CN---FORSYMMETRICCOLLOCATIONUSEDFORPARTICLEANDCM---FORASYMMETRICCOLLOCATIONUSEDFORCOLUMNN=7IW=1IS=2CALLCOLL(AS,BS,Q,XS,WS,19,N,IW,IS)NS=N+1WRITE(*,*)'*SYMMETRICSITUATION:*'WRITE(*,*)'*POLYNOMIALROOTS*'WRITE(*,*)(XS(I),I=1,NS)WRITE(*,*)WRITE(*,*)'*A-MATRIX*'DO20I=1,NS20WRITE(*,*)(AS(I,J),J=1,NS)WRITE(*,*)WRITE(*,*)'*B-MATRIX*'DO30I=1,NS30WRITE(*,*)(BS(I,J),J=1,NS)WRITE(*,*)WRITE(*,*)'*W-MATRIX*'WRITE(*,*)(WS(J),J=1,NS)CCALCULATINGTHEPARAMETERSOFTHEPROBLEM,WHICHWILLBEUSEDCFORTHEDIMENSIONLESSFORMOFANDDEFININGOFTHEPROBLEM.NEQ=2*NLRW=22+9*NEQ+NEQ**2LIW=20+NEQCINITIALCONDITIONSDO201I=1,NY(I)=1.D0Y(N+I)=0.D0201CONTINUEYN1=1.0D0YN2=0.D0T=0.D0DT=5.D-2ITOL=2RTOL=1.D-6DO203I=1,NEQATOL(I)=1.D-6203CONTINUEITASK=1ISTATE=1IOPT=0MF=22DO240IOUT=1,20TOUT=DT*DFLOAT(IOUT)CALLLSODE(FEX,NEQ,Y,T,TOUT,ITOL,RTOL,ATOL,ITASK,ISTATE,1IOPT,RWORK,LRW,IWORK,LIW,JEX,MF)OPEN(2,FILE='LW_S_ODE.OUT')WRITE(2,'(''Z:'',F8.4)')TWRITE(2,*)'R'WRITE(2,'(10(4X,D11.5))')(XS(I),I=1,N+1)WRITE(2,*)'T:'WRITE(2,'(10(4X,D11.5))')(Y(I),I=1,N),YN1WRITE(2,*)'C:'CDO205I=1,NWRITE(2,'(10(4X,D11.5))')(Y(N+I),I=1,N),YN2CWRITE(2,*)C205CONTINUEC220FORMAT(7HATT=,D12.4,6HY=,3D15.7)IF(ISTATE.LT.0)GOTO280240CONTINUEWRITE(3,260)IWORK(11),IWORK(12),IWORK(13)260FORMAT(/12HNO.STEPS=,I4,11HNO.F-S=,I4,11HNO.J-S=,I4)STOP280WRITE(3,290)ISTATE290FORMAT(///22HERRORHALT..ISTATE=,I3)STOPEND子程序：SUBROUTINEFEX(NEQ,T,Y,YP)CVARIABLESINTHEORDEROFY(1)TOY(M)FORTHECOLUMNCOLLOCATIONCPOINTSOF2TOM+1,Y(M+1)TOY(M+N)FORPARTICLECOLLOCATIONPOINTSCFROM1TONINCOLUMNCOLLOCATIONPOINT1,ANDY(M+N*(M+1)+1)TOCY(M+N*(M+2))FORPARTICLECOLLOCATIONPOINTSFROM1TONINCOLUMNCCOLLOCATIONPOINTM+2.IMPLICITREAL*8(A-H,O-Z)DIMENSIONY(19*19),YP(19*19),AS(19,19),BS(19,19)DOUBLEPRECISIONYN1,YN2COMMON/AB/N,AS,BSCOMMON/BC/YN1,YN2T0=0.D0TT0=0.D0DO990I=1,NT0=T0+AS(N+1,I)*Y(I)TT0=TT0+AS(N+1,I)*Y(N+I)990CONTINUEYN1=(0.92-T0)/(1+AS(N+1,N+1))YN2=-TT0/AS(N+1,N+1)DO99J=1,NT1=BS(J,N+1)*YN1T2=BS(J,N+1)*YN2DO991I=1,NT1=T1+BS(J,I)*Y(I)T2=T2+BS(J,I)*Y(I+N)991CONTINUEYP(J)=T1+0.2*(1-Y(N+J))*EXP(20-20/Y(J))YP(N+J)=T2+0.3*(1-Y(N+J))*EXP(20-20/Y(J))99CONTINUERETURNENDSUBROUTINEJEX(NEQ,T,Y,ML,MU,PD,NRPD)DOUBLEPRECISIONPD,T,YDIMENSIONY(NEQ),PD(NRPD,NEQ)RETURNEND输出结果：Z:.0500R....