如何解决这个悬而未决的难题?????

lonelyprince

我运行下面的程序出现了这个问题:
???  In an assignment  A(I) = B, the number of elements in B and
I must be the same.

Error in ==> cdpf_ABdif at 7
dy(3)=i/2*(Rabi1B*exp(i*dif*t)+Rabi2B*ExpZ)*y(4)-i*Rabi3*exp(i*dif*t)*y(2)+i/2*(Rabi1A+Rabi2A*ExpF)*y(1);

Error in ==> funfun\private\odearguments at 110
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ==> ode45 at 173
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...

Error in ==> CDPWnumerical_ABdif_rect at 62
[t,CDPWR]=ode45(@cdpf_ABdif,[t_begin:h_step:t_end],InitPara1,options);

我的主程序如下:
% Calulating population of Non-identical Two Tow-Level systems,
% with double Rectangular pulses, for CDPW and resonance
% By Xiang-yang Yu, XD Z, RZ Chen in SYSU
% 2008/04/10

clear all
close all
clc
global Rabi1A Rabi2A Rabi1B Rabi2B dif ExpZ ExpF

h=6.6260755e-34;                 % J s, 普朗克常数
hbar=h/2.0/pi;                      % J s，普朗克常数
mu2=9.73e-57;                   % (C m)^2, 电偶极矩平方   
R=0.616e-9;                     % m, 电偶极矩矩长

tau0=1.0e-15;                   % 飞秒的时间单位，s
Delta=0.05*1e+15;               % 失谐量，1/s   
Delta1=Delta*tau0;              % 无量纲的失谐量

%ccccccccccccccccccccccccccccccccccccccc
Tp1A0=50.0*1e-15;                 % s, Width of 1st Rectangular pulse with frequency A
Tp2A0=50.0*1e-15;                 % s, Width of 2nd Rectangular pulse with frequency A
Tp1B0=50.0*1e-15;                 % s, Width of 1st Rectangular pulse with frequency B
Tp2B0=50.0*1e-15;                 % s, Width of 2nd Rectangular pulse with frequency B
Tp1A=Tp1A0/tau0;                % fs, Width of 1st Rectangular pulse with frequency A
Tp2A=Tp2A0/tau0;                % fs, Width of 2nd Rectangular pulse with frequency A
Tp1B=Tp1B0/tau0;                % fs, Width of 1st Rectangular pulse with frequency B
Tp2B=Tp2B0/tau0;                % fs, Width of 2nd Rectangular pulse with frequency B
Area1A=pi;                      % 1st pulse area with frequency A
Area2A=0.5*pi;                  % 2nd pulse area with frequency A
Area1B=0.25*pi;                 % 1st pulse area with frequency B
Area2B=0.5*pi;                  % 2nd pulse area with frequency B
Rabi1A=Area1A/Tp1A;             % 1/fs, 1st Rabi frequency with frequency A
Rabi2A=Area2A/Tp2A;             % 1/fs, 2nd Rabi frequency with frequency A
Rabi1B=Area1B/Tp1B;             % 1/fs, 1st Rabi frequency with frequency B
Rabi2B=Area2B/Tp2B;             % 1/fs, 2nd Rabi frequency with frequency B
Rabi3=5*Rabi1A                  % 1/fs, controllable frequency induced by dipole-dipole interaction
dif=0.1                         % 1/fs, difference between frequency A and B

ExpZ=i;                         % exponential components introduced by pulses delay
ExpF=-i;

t_begin=0;                    
t_end=100;                       
h_step=0.1;
tspan=t_begin:h_step:t_end;     % time span

InitPara1=[1 0 0 0];            % The initial values of c1,c2,c3,c4
InitPara2=[0 1 0 0];
InitPara3=[0 0 1 0];
InitPara4=[0 0 0 1];

InitPara11=[0.5 0 0 0.5];       % The initial values of c1,c2,c3,c4
InitPara21=[0 0.5 0.5 0];
InitPara31=[0.5 0 0 0.5];
InitPara41=[0 0.5 0.5 0];

options=odeset('RelTol',1e-4,'AbsTol',1e-8);

% 相关双路径CDPW共振解
% while input=[1 0 0 0]
[t,CDPWR]=ode45(@cdpf_ABdif,[t_begin:h_step:t_end],InitPara1,options);

CDPWR_C1_1=(abs(CDPWR(:,1)).^2);
CDPWR_C2_1=(abs(CDPWR(:,2)).^2);
CDPWR_C3_1=(abs(CDPWR(:,3)).^2);
CDPWR_C4_1=(abs(CDPWR(:,4)).^2);

% while input=[0 1 0 0]
[t,CDPWR]=ode45('cdpf_ABdif',[t_begin:h_step:t_end],InitPara2,options);

CDPWR_C1_2=(abs(CDPWR(:,1)).^2);
CDPWR_C2_2=(abs(CDPWR(:,2)).^2);
CDPWR_C3_2=(abs(CDPWR(:,3)).^2);
CDPWR_C4_2=(abs(CDPWR(:,4)).^2);

% while input=[0 0 1 0]
[t,CDPWR]=ode45('cdpf_ABdif',[t_begin:h_step:t_end],InitPara3,options);

CDPWR_C1_3=(abs(CDPWR(:,1)).^2);
CDPWR_C2_3=(abs(CDPWR(:,2)).^2);
CDPWR_C3_3=(abs(CDPWR(:,3)).^2);
CDPWR_C4_3=(abs(CDPWR(:,4)).^2);

% while input=[0 0 0 1]
[t,CDPWR]=ode45('cdpf_ABdif',[t_begin:h_step:t_end],InitPara4,options);

CDPWR_C1_4=(abs(CDPWR(:,1)).^2);
CDPWR_C2_4=(abs(CDPWR(:,2)).^2);
CDPWR_C3_4=(abs(CDPWR(:,3)).^2);
CDPWR_C4_4=(abs(CDPWR(:,4)).^2);


hold on
subplot(2,2,1)
   plot(t,CDPWR_C1_1,'r', t,CDPWR_C1_2,'b', t,CDPWR_C1_3,'g', t,CDPWR_C1_4,'m');
   xlabel('Time (fs)');
   ylabel('|C_1|^2');
   axis([0 500 -0.1 1.1])
   set(gca,'Xtick',[0:50:500])
   legend('In=[1 0 0 0]','In=[0 1 0 0]','In=[0 0 1 0]','In=[0 0 0 1]');
   
subplot(2,2,2)
   plot(t,CDPWR_C2_1,'r', t,CDPWR_C2_2,'b', t,CDPWR_C2_3,'g', t,CDPWR_C2_4,'m');
   xlabel('Time (fs)');
   ylabel('|C_2|^2');
   axis([0 500 -0.1 1.1])
   set(gca,'Xtick',[0:50:500])
   
subplot(2,2,3)
   plot(t,CDPWR_C3_1,'r', t,CDPWR_C3_2,'b', t,CDPWR_C3_3,'g', t,CDPWR_C3_4,'m');
   xlabel('Time (fs)');
   ylabel('|C_3|^2');
   axis([0 500 -0.1 1.1])
   set(gca,'Xtick',[0:50:500])

subplot(2,2,4)
   plot(t,CDPWR_C4_1,'r', t,CDPWR_C4_2,'b', t,CDPWR_C4_3,'g', t,CDPWR_C4_4,'m');
   xlabel('Time (fs)');
   ylabel('|C_4|^2');
   axis([0 500 -0.1 1.1])
   set(gca,'Xtick',[0:50:500])

待求解的微分方程为:
%相关双路径运动方程（微分方程组）
%with dipole-dipole interaction

function dy = cdpf_ABdif(t,y)
global Rabi1A Rabi2A Rabi1B Rabi2B Rabi3 ExpZ ExpF dif
dy=zeros(4,1)         
dy(4)=i/2*(Rabi1B+Rabi2B*ExpF)*y(3)+i/2*(Rabi1A+Rabi2A*ExpF*exp(i*dif*t))*y(2);
dy(3)=i/2*(Rabi1B*exp(i*dif*t)+Rabi2B*ExpZ)*y(4)-i*Rabi3*exp(i*dif*t)*y(2)+i/2*(Rabi1A+Rabi2A*ExpF)*y(1);
dy(2)=i/2*(Rabi1A+Rabi2A*ExpZ*exp(-i*dif*t))*y(4)-i*Rabi3*exp(-i*dif*t)*y(3)+i/2*(Rabi1B*exp(-i*dif*t)+Rabi2B*ExpF)*y(1);
dy(1)=i/2*(Rabi1A+Rabi2A*ExpZ*exp(-i*dif*t))*y(3)+i/2*(Rabi1B*exp(i*dif*t)+Rabi2B*ExpZ)*y(2);

我是新手, 许多细节的问题都不明白, 还望各位高手帮忙, 先谢了!!!!

ChaChing

我设断点看了下!
cdpf_ABdif函数中dy(3)式中Rabi3为空矩阵! dy(3)=[]有问题!
回头找发现LZ主程序与cdpf_ABdif函数间的global不同, 具体应如何, LZ应最清楚!

lonelyprince

问题已解决, 正如楼主所言! Thanks a lot!