function my_ode()
%parameters used inthe simulation.
pi= 3.143; lambda_1 = 1536 ; lambda_2 = 1538; lambda_s = 1534; lambda_0 = 1535;
lambda_c = (lambda_1 + lambda_2)/2;
w_0 = (2*pi)/lambda_0;
w_1 = (2*pi)/lambda_1;
w_2 = (2*pi)/lambda_2;
w_s = (2*pi)/lambda_s;
w_c = (2*pi)/lambda_c;
B_3 =2.881*10^-2;B_4 =-1.852*10^-4;P_1 = 1;P_2 =1 ;
gamma =122.7*10^-3;
B_2 = B_3*(w_c - w_0)^2 + (B_4/2)*(w_c - w_0);
delta_w_p = (w_2 - w_1)/2;
delta_w_s = w_s - w_c;
w_sb1 = 2*delta_w_p - delta_w_s;
%w_sb2 = 2*delta_w_p - delta_w_s;
delta_w_sb = w_c - w_sb1;
r = 2*gamma*sqrt(P_1*P_2);
A_sb10 = 0;
A_s0 = 1;
A_10 = 0;
A_sb20 = 0;
% function to compute derivative of
%A = [A_s, A_1, A_sb1, A_sb2];
function dAdz = odefcn(~,A)
dA_sb1dz = (1i/24)*A_sb1*(12*B_2*(delta_w_sb^2 - delta_w_p^2)...
- 4*B_3*(delta_w_s^3-delta_w_p^3)+ B_4*(delta_w_sb^4 - delta_w_p^4))+1i*gamma...
*P_1*(A_sb1 + conj(A_s))+ 1i*r*(conj(A_sb1) + A_s);
dA_sdz = (1i/24)*A_s(12*B_2*(delta_w_s^2 - delta_w_p^2)...
- 4*B_3*(delta_w_s^3-delta_w_p^3)+ B_4*(delta_w_s^4 - delta_w_p^4))+ 1i*gamma...
*P_1*(conj(A_sb1) + A_s)+ 1i*r*(conj(A_1) + A_sb2);
dA_1dz = (1i/24)*A_1(12*B_2*(delta_w_s^2 - delta_w_p^2)...
+ 4*B_3*(delta_w_s^3-delta_w_p^3)+ B_4*(delta_w_s^4 - delta_w_p^4))+1i*gamma...
*P_2*(conj(A_sb2) + A_1)+ 1i*r*(conj(A_s) + A_sb1);
dA_sb2dz = (1i/24)*A_sb2(12*B_2*(delta_w_sb^2 - delta_w_p^2)...
+ 4*B_3*(delta_w_sb^3-delta_w_p^3)+ B_4*(delta_w_sb^4 - delta_w_p^4))+1i*gamma...
*P_2*(conj(A_1) + A_sb2)+ 1i*r*(conj(A_sb1) + A_s);
dAdz = A;
dAdz(:) = [dA_sb1dz dA_sdz dA_1dz dA_sb2dz];
end
A0 = [A_s0 A_10 A_sb10 A_sb20];
[z,A] = ode45(@odefcn, [0 1], A0);
A_sb1 = A(:,1); A_1 = A(:,3);
A_s = A(:,2); A_sb2 = A(:,4);
end

3 commentaires

Select the whole code and then press the code button {}
kumbirayi
kumbirayi le 18 Oct 2018
thank you
input datas??

Connectez-vous pour commenter.

 Réponse acceptée

Torsten
Torsten le 18 Oct 2018
Please check whether this is what you want:
function main
A_sb10 = 0;
A_s0 = 1;
A_10 = 0;
A_sb20 = 0;
A0 = [A_s0 A_10 A_sb10 A_sb20];
[z,A] = ode45(@(t,y)odefcn(t,y(1),y(2),y(3),y(4)), [0 1], A0);
A_sb1 = A(:,1); A_1 = A(:,3);
A_s = A(:,2); A_sb2 = A(:,4);
plot(z,real(A(:,4)))
end
function dAdz = odefcn(t ,A_s, A_1, A_sb1, A_sb2)
%parameters used inthe simulation.
pi= 3.143; lambda_1 = 1536 ; lambda_2 = 1538; lambda_s = 1534; lambda_0 = 1535;
lambda_c = (lambda_1 + lambda_2)/2;
w_0 = (2*pi)/lambda_0;
w_1 = (2*pi)/lambda_1;
w_2 = (2*pi)/lambda_2;
w_s = (2*pi)/lambda_s;
w_c = (2*pi)/lambda_c;
B_3 =2.881*10^-2;B_4 =-1.852*10^-4;P_1 = 1;P_2 =1 ;
gamma =122.7*10^-3;
B_2 = B_3*(w_c - w_0)^2 + (B_4/2)*(w_c - w_0);
delta_w_p = (w_2 - w_1)/2;
delta_w_s = w_s - w_c;
w_sb1 = 2*delta_w_p - delta_w_s;
%w_sb2 = 2*delta_w_p - delta_w_s;
delta_w_sb = w_c - w_sb1;
r = 2*gamma*sqrt(P_1*P_2);
dA_sb1dz = (1i/24)*A_sb1*(12*B_2*(delta_w_sb^2 - delta_w_p^2)...
- 4*B_3*(delta_w_s^3-delta_w_p^3)+ B_4*(delta_w_sb^4 - delta_w_p^4))+1i*gamma...
*P_1*(A_sb1 + conj(A_s))+ 1i*r*(conj(A_sb1) + A_s);
dA_sdz = (1i/24)*A_s*(12*B_2*(delta_w_s^2 - delta_w_p^2)...
- 4*B_3*(delta_w_s^3-delta_w_p^3)+ B_4*(delta_w_s^4 - delta_w_p^4))+ 1i*gamma...
*P_1*(conj(A_sb1) + A_s)+ 1i*r*(conj(A_1) + A_sb2);
dA_1dz = (1i/24)*A_1*(12*B_2*(delta_w_s^2 - delta_w_p^2)...
+ 4*B_3*(delta_w_s^3-delta_w_p^3)+ B_4*(delta_w_s^4 - delta_w_p^4))+1i*gamma...
*P_2*(conj(A_sb2) + A_1)+ 1i*r*(conj(A_s) + A_sb1);
dA_sb2dz = (1i/24)*A_sb2*(12*B_2*(delta_w_sb^2 - delta_w_p^2)...
+ 4*B_3*(delta_w_sb^3-delta_w_p^3)+ B_4*(delta_w_sb^4 - delta_w_p^4))+1i*gamma...
*P_2*(conj(A_1) + A_sb2)+ 1i*r*(conj(A_sb1) + A_s);
dAdz = [dA_sb1dz; dA_sdz; dA_1dz; dA_sb2dz];
end

1 commentaire

kumbirayi
kumbirayi le 19 Oct 2018
Modifié(e) : kumbirayi le 19 Oct 2018
it does make sense, i can work with this code. Thank you!

Connectez-vous pour commenter.

Plus de réponses (0)

Catégories

En savoir plus sur Numerical Integration and Differential Equations dans Centre d'aide et File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by