my code is entering infinite loop when i am trying to solve differential equation using ode45 command , can some one help me fix this
3 views (last 30 days)
Show older comments
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);
3 Comments
Walter Roberson
on 10 Oct 2018
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.
Accepted Answer
Stephan
on 10 Oct 2018
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
2 Comments
Stephan
on 11 Oct 2018
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
More Answers (0)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!