my code is entering infinite loop when i am trying to solve differential equation using ode45 command , can some one help me fix this

 function dE_dz =asim(z,E1,~,~)

 c=3*1e8;
L=800*10^-9
pi=3.1415926535;
Dk=20.94;
w1=(2*pi*c)/L;
%w1=fs*[-N/2:N/2-1];
deff=2.0657e-12;
Lc=2*pi/Dk;
Lnl=50*Lc;
n2=5*10^-16;
I0=10^-9;

 dE_dz = (-1i*conj(E1).*E1.^2.*w1*deff)/(c*n2*Dk)+ 1i*2*pi*(n2*I0)*(Lnl/L)*((abs(E1)).^2).*E1 +...
    + 1i*4*pi*(n2*I0)*(Lnl/L).*((-w1*deff.*E1.^2*exp(-1i*Dk*z))/(c*n2*Dk)).^2.*E1; 
end

 clc; clear all; close all; clf;

 Pi=1.98;                               %input('Enter the value of input power in mW  ')
t=120*10^-15                           %input('Enter the value of input pulse width in seconds  ')
tau=-300*10^-15:0.1*10^-15:300*10^-15;                         %input('Enter time period with upper(U), lower(L) and interval between upper and  lower interval(I) in this format L:I:U')
dt=1e-15/t
fs=1/dt;
pi=3.1415926535;
Ao=sqrt(Pi)                            % amplitude of an input gaussian pulse
h=2000                                 % small step as defined by  split step algorithm
                                       %for ii=0.1:0.1:(s/10)              %1*10^-7:1*10^-7:s %different fiber lengths
E1=Ao*exp(-(tau/t).^2); % generation of an gaussian input pulse
figure(1)
plot(tau,abs(E1));                         % graph of input pulse
title('Input Pulse'); xlabel('Time in ps'); ylabel('Amplitude');
grid on;
hold on;

 [z,E1]=ode45(@asim,[0 1.0],0.6325)
%q = ode45(@(z,E1)asim(z,E1,~,~),[0 1],E1_0);
% time=
% % % tau1=L*(-300*10^-15:300*10^-15);
plot(z,E1);

[madhan ravi]

upload asim function and format your code correctly to get an answer

[asim asrar]

i have uploaded the asim function prior to the code , i will update my code but i am not getting where i am doing mistake , please help me with that

[Walter Roberson]

Why do you think you have an infinite loop? Perhaps your equations are just stiff or oscillate enough that they have to be integrated with a small step size.

Answer 1

[Stephan]

Hi,

You make your code inefficient by calculating constant values again and again in every iteration. Try this:

syms E1 z
c=3*1e8;
L=800*10^-9;
pi=3.1415926535;
Dk=20.94;
w1=(2*pi*c)/L;
deff=2.0657e-12;
Lc=2*pi/Dk;
Lnl=50*Lc;
n2=5*10^-16;
I0=10^-9;
dE_dz = (-1i*conj(E1).*E1.^2.*w1*deff)/(c*n2*Dk)+ 1i*2*pi*(n2*I0)*(Lnl/L)*((abs(E1)).^2).*E1 +...
    + 1i*4*pi*(n2*I0)*(Lnl/L).*((-w1*deff.*E1.^2*exp(-1i*Dk*z))/(c*n2*Dk)).^2.*E1; 
matlabFunction(dE_dz,'File','asim_fast')

This code creates a file asim_fast.m containing this code - which does the same calculation much faster:

function dE_dz = asim_fast(E1,z)
%ASIM_FAST
%    DE_DZ = ASIM_FAST(E1,Z)

%    This function was generated by the Symbolic Math Toolbox version 8.2.
%    10-Oct-2018 20:02:27

t2 = E1.^2;
t3 = t2.^2;
t4 = abs(E1);
dE_dz = E1.*t4.^2.*5.891597660294395e-17i-t2.*conj(E1).*1.549567321951994e9i+E1.*t3.*exp(z.*-4.188e1i).*2.829332414080319e2i;

Change your ode call to:

Pi=1.98;
t=120*10^-15;
tau=-300*10^-15:0.1*10^-15:300*10^-15;
dt=1e-15/t;
fs=1/dt; 
Ao=sqrt(Pi);
h=2000;
E1=Ao*exp(-(tau/t).^2);
figure(1)
subplot(2,1,1)
plot(tau,abs(E1));
title('Input Pulse');
xlabel('Time in ps');
ylabel('Amplitude'); 
grid on; 
[z,E1]=ode45(@asim_fast,[0 1.0],0.6325);
subplot(2,1,2)
plot(z,E1)

This will significant increase speed.

I used subplot to better see the both diagrams. Since the result is changing mainly in the imaginary part perhaps another kind of plot (e.g. z over abs(E1) would make more sense - but i have no insight to your problem.

Best regards

Stephan

[asim asrar]

Dear Stephan, Thenks for the assistance. but when i am trying to manipulate the code according to you its giving. Error using sym/matlabFunction>writeMATLAB (line 382) Could not create file asim_fast.m: Permission denied

Error in sym/matlabFunction (line 151) g = writeMATLAB(funs,file,varnames,outputs,body, opts.optimize); what do i need to do further

[Stephan]

You could copy the content i posted and save the file by hand. Or change the line where the filename is specified to a path you are allowed to write in. Alternativly use in one file:

Pi=1.98;
t=120*10^-15;
tau=-300*10^-15:0.1*10^-15:300*10^-15;
dt=1e-15/t;
fs=1/dt; 
Ao=sqrt(Pi);
h=2000;
E1=Ao*exp(-(tau/t).^2);
figure(1)
subplot(2,1,1)
plot(tau,abs(E1));
title('Input Pulse');
xlabel('Time in ps');
ylabel('Amplitude'); 
grid on; 
[z,E1]=ode45(@asim_fast,[0 1.0],0.6325);
subplot(2,1,2)
plot(z,E1)

function dE_dz = asim_fast(E1,z)
%ASIM_FAST
%    DE_DZ = ASIM_FAST(E1,Z)
%    This function was generated by the Symbolic Math Toolbox version 8.2.
%    10-Oct-2018 20:02:27
t2 = E1.^2;
t3 = t2.^2;
t4 = abs(E1);
dE_dz = E1.*t4.^2.*5.891597660294395e-17i-t2.*conj(E1).*1.549567321951994e9i+E1.*t3.*exp(z.*-4.188e1i).*2.829332414080319e2i;
end