[交流] 【求助】一个脉冲热效应的matlab程序不能运行

本人编写了一个脉冲热效应的matlab 程序，显示结果为边界值，也就是主程序不能运行，即 if=1后面出现了问题，恳求有关高手给予指点。程序如下 lear clc % ************************************* a=1.5;b=5;%晶体参数 h=0.1;ht=5e-4;%计算步长 tol=0.001;maxl=3000;%***误差 kr=0.014;%热导率 absorb=0.91;%晶体吸收系数 w0=0.32;%泵光半径 Q=15;%泵浦功率 dndt=7.3e-6;%dn/dt n0=1.82;%折射率 o=7.5e-6;%热膨胀系数 L=30;%腔长 R=50;%输出镜半径 TOUT=0.1;%输入镜透过率 ls=946e-6;lp=808e-6;lf=1062e-6;%激光波长，泵浦波长，平均荧光波长 oabs=7.9e-24;oem=4e-24;%受激吸收截面，受激发射截面 t=230e-6;%上能级寿命 tl=0.025;%脉冲持续时间 rzxs=4e-4;   % when adopt inter must follow  inter<=f, inter=f maybe the best c=3e11;%光速 hc=6.62e-34;%普朗克常数 n=fix(a/h)+1;m=fix(b/h)+1;nt=fix(tl/ht)+1; % importent nt=fix(tl/ht)+1; otherwise the number will decrease inter%计算取用格子数 q=zeros(2*n-1,m); % ******************************** h0=0.006;%对流换热系数 %********************************* % 边界条件 Es=hc*c/ls;Ep=hc*c/lp;%光子能量 K=0.006;%光光转化效率 fa=0.074;%玻耳兹曼因子 lr=0.015;%腔内损耗 ql=4.22; %晶体密度 cl=0.59;%晶体比热 ws=((ls/pi)^2*L*(R-L))^(0.25);%激光半径公式 z0=n0*w0^2*pi/lp; for j=1:m     wp(j)=w0*(1+j*h/z0)^(0.5);%泵浦光斑半径公式     pp(j)=Q*exp(-absorb*j*h);%泵浦光传播功率      ps=Q*K/TOUT*exp(0.5);%激光传播功率 end %程序主体 U=ones(2*n-1,m,nt); U0=16; U(1,1:m,1:nt)=16;U(2*n-1,1:m,1:nt)=16;U(1:2*n-1,1,1:nt)=16;U(1:2*n-1,m,1:nt)=16;U(1:2*n-1,1:m,1)=16; ******************************************************           for T=2:nt              for i=2:2*n-2                  for j=1:m                       q(i,j)=2*Q*(1-(lp/lf+(lp/ls-lp/lf)*...           (ps*exp(-2*(wp(j)/2)^2/ws^2)*oem/Es)/(1+(t*ps*2/(pi*ws^2)*exp(-2*(wp(j)/2)^2/ws^2)*oem/Es))*...      (t-fa/(pp(j)*2/(pi*wp(j)^2)*exp(-2*(wp(j)/2)^2/ wp(j)^2)*oabs/Ep))))*absorb/(3.14*w0^2*(1-exp(-absorb*b)))*exp(-2*((i-n)*h)^2/w0^2)*exp(-absorb*j*h);                  if j==1                            U(i,j,T)=U(i,j,T-1)+kr*ht/ql/cl*((U(i+1,j,T-1)-2*U(i,j,T-1)+U(i-1,j,T-1))/h^2+(U(i,j,T-1)-U(i-1,j,T-1))/(i-n+1)/h^2+(-2*U(i,j,T-1)-2*h*rzxs*(U(i,j,T-1)-U0)+2*U(i,j+1,T-1)))+q(i,j)*ht/ql/cl                                                     elseif j==m                                                            U(i,j,T)=U(i,j,T-1)+kr*ht/ql/cl*((U(i+1,j,T-1)-2*U(i,j,T-1)+U(i-1,j,T-1))/h^2+(-U(i-1,j,T-1))/(i-n+1)/h^2+(U(i,j+1,T-1)-2*U(i,j,T-1)+2*h*rzxs*(U(i,j,T-1)-U0)+2*U(i,j-1,T-1)))+q(i,j)*ht/ql/cl                          else                              U(i,j,T)=U(i,j,T-1)+kr*ht/ql/cl*((U(i+1,j,T-1)-2*U(i,j,T-1)+U(i-1,j,T-1))/h^2+(U(i,j,T-1)-U(i-1,j,T-1))/(i-n+1)/h^2+(U(i,j+1,T-1)-2*U(i,j,T-1)+U(i,j-1,T-1)))+q(i,j)*ht/ql/cl                                                   end                  end                end          end % ******************************************** output=U(n+1,5,; output=output(; plot(output) [ Last edited by nono2009 on 2009-12-24 at 22:02 ]

回复此楼