% Fig_4_5.m % classical oscillator using the RK4 algorithm to see time evolution. % The different physical condtions can be achived by adjusting the % parameters. Up to 3rd-order non-linearity is included. % Linear term is linear_coeff*(omegax)^2 and 3rd-order term is % non_linear_3*(omegax)^2. %plotting parameters + fontsizes clear all; clf; FS = 28; %label fontsize 18 FSN = 28; %number fontsize 16 LW = 2; %linewidth 2 % Change default axes fonts. set(0,'DefaultAxesFontName', 'Times New Roman') set(0,'DefaultAxesFontSize', FSN) % Change default text fonts. set(0,'DefaultTextFontname', 'Times New Roman') set(0,'DefaultTextFontSize', FSN) %%%%%%%%%%%%%%%% input parameters data %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% omegax=pi*1e+06; % unit:s-1, frequency of the oscillator. 3.2698e+06 Hz domegax=omegax*0.10; % frequency shift (0, 0.01, 0.1, 1.0) omegax1=(pi*1e+06)+domegax; % unit:s-1, frequency of the periodic force % applied to the oscillator. 3.2691e+06 gamma = omegax / 200;% unit:s-1, oscillator damping term (200, 500) linear_coeff =-1.10; % linear simple harmonic term, +ve (1) for SHO % and -ve (-1, -1.5, -2.5) for Duffing oscillator non_linear_3 =1.0; % coefficient of 3rd order nonlinear term. (0, 1, 10) %Example of SHO use /200, 1.0, 0.0, 50e12 %Example frequency doubling use /200, -2.5, 1.0, 50e12 total_time=0.01; % unit:s, total time for which system is observed. 0.1 fraction1= 0.01; % time increment of integration = fraction1 / omegax %it is a small fraction of the oscillatory period %set to smaller value for smaller time steps. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% enter initial values %%%%%%%%%%%%%%% y1= 0; % initial VELOCITY of oscillator at t = 0.0 y2= 0; % initial POSITION of oscillator at t = 0.0 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% deltat1=(fraction1*2*pi/omegax);% unit:s, time increment of integration, % set to 1% of oscillation period deltat2=(fraction1*2*pi/omegax1);% unit:s, time increment of integration, % set to 1% of oscillation period %smaller of the time increments chosen for integretion. arr = sort([deltat1 deltat2]); deltat = arr(1); extent=round(total_time/deltat) % total number of time steps determined % from time increment. ss = 1; if gamma>0 ss = round((50/gamma)/deltat) % remove initial transient decay % from position array end F0(1:extent) = 50e12; % amplitude of periodic force applied to oscillator % in absence of noise. 1e11 - 1e14 % F0(ss+round(0.005*extent)-1:extent)=-6335e12;%-6335e12 % F0(ss+round(0.005*extent))=+6335e12;%-6335e12 % F0(ss+round(0.005*extent):extent)=-6340e12;%-6335e12 % F0(ss+round(0.005*extent)+50:extent)=0.0; %%%%%%%%%%%%% initialising storage matrices%%%%%%%%%%%%%%%%%%%%%%%%%% Y(extent,2)=0;% array with 4 coloumns, column 1 : velocity 1, % column 2 : position 1, column 3 : velocity 2, column 4 : position 2. %T(1:extent)=0; % stores time values at every integration point. T=ones(1,extent); %initialize time array T=deltat*double(find(T));%initialize time values at every integration point lcoeff=linear_coeff*(omegax^2);% redefine coeff_of_linear nlcoeff =non_linear_3*(omegax^2);% redefine 3rd-order nonlinear term Force(1:extent)=0; % stores force values at every integration point. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%% integration loop %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for index=1:extent %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% RK4 algorithm %%%%%%%%%%%%%%%%%%%%%%%%%% Force(index) = (F0(index))*cos(omegax1*T(index)) ; % update force array temp1=y1;%velocity temp2=y2;%position % 1st step F = (F0(index))*cos(omegax1*T(index)) ; d1 = -lcoeff*y2 - nlcoeff*(y2^3) - gamma*y1 + F;%cos(omegax1*t); d2 = y1; ak1=d1*deltat; bk1=d2*deltat; y1=temp1+0.5*ak1; y2=temp2+0.5*bk1; %2nd step F = (F0(index))*cos(omegax1*(T(index)+deltat/2)); d1 = -lcoeff*y2 - nlcoeff*(y2^3) - gamma*y1+F;%cos(omegax1*t); d2 = y1; ak2=d1*deltat; bk2=d2*deltat; y1=temp1+0.5*ak2; y2=temp2+0.5*bk2; %3rd step F = (F0(index))*cos(omegax1*(T(index)+deltat/2)) ; d1 = -lcoeff*y2 - nlcoeff*(y2^3) - gamma*y1+F;%cos(omegax1*t); d2 = y1; ak3=d1*deltat; bk3=d2*deltat; y1=temp1+1*ak3; y2=temp2+1*bk3; %4th step F = (F0(index))*cos(omegax1*(T(index)+deltat)) ; d1 = -lcoeff*y2 - nlcoeff*(y2^3) - gamma*y1+F;%cos(omegax1*t); d2 = y1; ak4=d1*deltat; bk4=d2*deltat; y1 = temp1+(ak1+2.0*ak2+2.0*ak3+ak4)/6.0; y2 = temp2+(bk1+2.0*bk2+2.0*bk3+bk4)/6.0; % store position and velocity values. Y(index,1)=y1; Y(index,2)=y2; end save Y.mat Y save T.mat T %Plotting the columns of the returned array Y versus T shows the solution ttl=(['lin=',num2str(linear_coeff),', nonlin=',... num2str(non_linear_3),', f_0 =',num2str(omegax/1e6, '%7.2f')... ,' MHz, f_F =',num2str(omegax1/1e6, '%7.2f')... ' MHz, \gamma =',num2str(gamma/1e6, '%7.4f'),' MHz']); figure(1) plot(T(ss:ss+round(0.001*extent))*1e6/(2*pi),... (Y(ss:ss+round(0.001*extent),2)),'k','LineWidth',LW); xlabel('Time, \it{t}\rm (\mus)','FontSize',FS); ylabel('Position, \it{x}\rm','FontSize',FS); axis([T(ss)*1e6/(2*pi),T(ss+round(0.001*extent))*1e6/(2*pi),... 1.1*min(Y(:,2)),1.1*max(Y(:,2))]); %title(ttl,'FontSize',FS); set(gca,'LineWidth',LW) set(gca,'TickLength',[0.02 0.025]) figure(2) plot(Y(ss:end,2)/1000,(Y(ss:end,1)/1000),'k','LineWidth',LW); xlabel('Position, \it{x}\rm(\it{t}\rm)','FontSize',FS); ylabel('Velocity, d\it{x}\rm/d\it{t}\rm','FontSize',FS); %title(ttl,'FontSize',FS); set(gca,'LineWidth',LW) set(gca,'TickLength',[0.02 0.025]) av_amp_osc1=sum(abs(Y(:,2)))/extent max_amp_osc1=max(abs(Y(:,2))) %calculate total energy in system Total_energy=0.5*(Y(:,1).^2+linear_coeff*(omegax^2)*Y(:,2).^2) +... (non_linear_3/4)*(omegax^2)*Y(:,2).^4 ; axis ([-4e-3 4e-3 -3e4 3e4]); max_amp = 1.1*max(abs(Y(:,2))); potential_position = -max_amp:0.001*max_amp:max_amp; Potential = 0.5*linear_coeff*(omegax^2)*potential_position.^2 + ... (non_linear_3/4)*(omegax^2)*potential_position.^4 ; xmin=min(potential_position); xmax=max(potential_position); ymin=min(Potential); ymax=max(Potential); figure(3) plot(potential_position,Potential,'k','LineWidth',LW); xlabel('Position, \it{x}\rm','FontSize',FS); ylabel('Potential, \itV\rm(\it{x}\rm)','FontSize',FS); axis([xmin/2,xmax/2,ymin,0.03*ymax]); % title(['Potential, lin=',num2str(linear_coeff),', nonlin=',... % num2str(non_linear_3),', \omega =',num2str(omegax, '%7.2e')]); set(gca,'LineWidth',LW) set(gca,'TickLength',[0.02 0.025])