Impostare equazioni su ODE45

[nandotp]

Mi sta venendo una confusione pazzesca

[dimaios]

Non c'è alcuna confusione.

Alcuni libri riportano la seconda equazione con il segno meno nel primo termine ed altri con il segno più.
Bisogna capirne il motivo, al limite puoi anche chiedere al docente facendo vedere vari riferimenti.

Guardando le dimostrazioni su come derivano il risultato puoi capire dov'è il problema ; ovviamente è un lavoro piuttosto lungo ed impegnativo.

[nandotp]

Comunque secondo me i segni sono quelle che ci sono nelle slide perché se lei da una occhiata veloce nel file raman scattering.pdf (da pag 14 a 15) capisce il motivo. Ho visto in altri libri che sono entrambi negativi i segni.
Breve spiegazione:
Il primo addendo nel secondo membro esprime la quantità di atomi che decadono per emissione stimolata mentre il secondo esprime la quantità di atomi che si eccitano per assorbimento stimolato.

[dimaios]

Ho visto in altri libri che sono entrambi negativi i segni.

Ed in altri sono invece discordi
E' inutile girarci intorno, questo è un problema da chiarire.

Comunque per non perdere tempo puoi simulare utlizzando le equazioni semplificate per vedere se il risultato è coerente.

[nandotp]

Facendo la simulazione con l'equazioni semplificate ottengo i seguenti grafici:
Che sembrano più coerenti con gli andamenti che dovrei ottenere ma se si guarda il grafico dell'uscita nel post 5 con Po=1 Watt e Ps0=1mW non è proprio uguale, Io ho plottato il guadagno usando semilogy. Ancora sono in alto mare perché tutti i grafici divergono. Il codice matlab è il seguente:

Codice: Seleziona tutto

function dy=Raman2(z,y)
%Co-propagating system equations
dy=zeros(2,1);


%y(1)==Is Intensity input Signal
%y(2)==Ip Intensity of the pump
%Is_dot=dy(1)=Is'=dIs/dz
%Ip_dot=dy(2)=Ip'=dIp/dz
Aeff=2.5*10^-9; %m^2
g_r=10.^-13; %coeff. Raman [m/W]
Lambda_p=1550*10^(-9); %Wavelength of pump [nm]
Lambda_s=1450*10^(-9);%Wavelength signal [nm]
fp=3e8/Lambda_p; % Pump Frequency [Hz]
 fs=3e8/Lambda_s; %Signal Frequency  [Hz]
 f_ratio=fp/fs;


alpha_s=0.046; %attenuation [km.^-1]
alpha_p=alpha_s;
dy(1)=g_r.*y(1).*y(2);
dy(2)=-f_ratio.*g_r.*y(1).*y(2);
 
end

Codice: Seleziona tutto

%SOLUZIONE NUMERICA
 
clear all
close all

P_p0=1; %Potenza pompa iniziale [W]
P_s0=0.001;%Potenza segnale iniziale [W] o [0 dBm]

% P_s0_dBm=10*log10(P_s0/0.001);
% P_p0_dBm=10*log10(P_p0/0.001);

span=[0:0.1:5000]; %lunghezza tratta [m]
Aeff=(80*10.^-6).^2; %area efficace micro_m^2

 Is_0=P_s0/Aeff;
 Ip_0=P_p0/Aeff;

y1=[Is_0];
y2=[Ip_0];
y0=[y1;y2]; %condizioni iniziali P_s0=Ps(Z=0)=1mW ; P_p0=1W potenza pompa iniziale

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

[z,y]=ode45('Raman2',span,y0,options);
figure
plot(z,y(:,1),'r')
legend('Is')
xlabel('Lunghezza [m]')
ylabel('Densità di Potenza segnale [W/(m^2)]')
title('Andamento Intensità segnale')
grid on



figure
plot(z,y(:,2),'b')
legend('Ip')
xlabel('Lunghezza [m]')
ylabel('Densità di Potenza Pompa [W/(m^2)]')
title('Andamento Intensità Pompa')
grid on

G=10*log10(y(:,1)/Is_0);

figure
semilogy(z,G*10.^3,'m')
legend('Guadagno')
xlabel('Lunghezza [m]')
ylabel('Guadagno ')
title('Andamento Guadagno Raman')
grid on


[dimaios]

Ti suggerivo di graficare le soluzioni delle equazioni differenziali in forma chiusa e semplificata già risolte nelle dispense.

[nandotp]

Ho scritto il codice relativo alla forma chiusa ma non torna niente. Ecco i grafici e il codice (spero di non aver sbagliato niente).

Codice: Seleziona tutto

clear all
close all

P_p0=1; %Potenza pompa iniziale [W]
P_s0=0.001;%Potenza segnale iniziale [W] o [0 dBm]
span=[0:1:5000]; %lunghezza tratta [m]
Aeff=(80*10.^-6).^2; %area efficace micro_m^2
 P_p=[1 2 3]; %vettore potenze pompa [W]
Ip=P_p./Aeff; %vettore intensità di pompa con Pp0=1W,Pp1=2W e Pp2=3W
g_r=10.^(-13); %Coefficiente Raman [m/W]
% P_s0_dBm=10*log10(P_s0/0.001);
% P_p0_dBm=10*log10(P_p0/0.001);

 Is_0=P_s0/Aeff;%intensità segnale [W/m^2]

 

 for i=1:length(span)
     for j=1:length(Ip)
   Is(i,j)=Is_0*exp(g_r*Ip(1,j).*span(i));
     end
 end
 
 
 figure
plot(span,Is(:,1),span,Is(:,2),span,Is(:,3))
legend('Intensità segnale 1 con Pp0=1W ','Intensità segnale 2 con Pp0=2W ','Intensità segnale 3 con Pp0=3W')
xlabel('Lunghezza [m]')
%  ylim([0 30])
ylabel('Ampiezza [W/(m^2)]')
title('Andamento intensità segnale in fibra')
grid on

Ps1=Aeff.*Is(:,1);
Ps2=Aeff.*Is(:,2);
Ps3=Aeff.*Is(:,3);
Ps=[Ps1 Ps2 Ps3];
Ps_z_dB=10*log10(Ps*10^3);

%P[dBm]=10*log10(P[mW]/1mW);

figure
semilogy(span,Ps_z_dB(:,1),span,Ps_z_dB(:,2),span,Ps_z_dB(:,3))
legend('Potenza segnale 1 con Pp0=1W ','Potenza segnale 2 con Pp0=2W ','Potenza segnale 3 con Pp0=3W')
xlabel('Lunghezza [m]')
%  ylim([0 30])
ylabel('Ampiezza [dBm]')
title('Andamento della potenza del segnale in fibra')
grid on


[dimaios]

L'area efficace è pari a

Codice: Seleziona tutto

Aeff=(80*10.^-6).^2; %area efficace micro_m^2


Tre errori in una riga di programma ; uno nella sostanza uno di programmazione e uno nel commento.
Qui si concludono i miei interventi in questo thread.

[nandotp]

Buon pomeriggio, mi scuso per l'errore che ho commesso, dovevo scrivere per correttezza . Penso di aver risolto il problema, in pratica avevo scambiato la lunghezza d'onda della pompa con quella del segnale , un altro errore è che l'intensità della pompa deve essere maggiore del
segnale e non il viceversa, dopodichè ho modificato in scala lineare dividendolo per 1000.
Posto il seguente risultato:
In questo grafico vi è il confronto tra la soluzione con ode45 e quella esatta. Vedendo le slide mi sembrano quasi uguali.
Allego anche il codice matlab:

Codice: Seleziona tutto

%SOLUZIONE NUMERICA
 
clear all
close all
span=[0:1:5000]; %lunghezza tratta [m]

Aeff=80*10.^-12; %area efficace um^2


P_p0=1; %Potenza pompa iniziale [1W]
P_p1=2;  %Potenza pompa iniziale [2W]
P_p2=3;  %Potenza pompa iniziale [3W]

Is_0=0.001/Aeff;%1mW/m^2
Is_1=0.002/Aeff;%2mW/m^2
Is_2=0.003/Aeff;%3mW/m^2

Ip_0=P_p0/Aeff;
Ip_1=P_p1/Aeff;
Ip_2=P_p2/Aeff;

y0=[Is_0;Ip_0]; %condizioni iniziali1
y1=[Is_1;Ip_1];%condizioni iniziali2
y2=[Is_2;Ip_2];%condizioni iniziali3

options= odeset('RelTol',1e-12,'AbsTol',1e-9);
save options

[z,y_0]=ode45('Raman2',span,y0,options);
[z,y_1]=ode45('Raman2',span,y1,options);
[z,y_2]=ode45('Raman2',span,y2,options);
g_r=10.^(-13); %coefficiente Raman [m/W]


figure
plot(z,y_0(:,1),'r')
legend('Is0')
xlabel('Lunghezza [m]')
ylabel('Densità di Potenza segnale 1[W/(m^2)]')
title('Andamento Intensità segnale 1')
grid on

figure
plot(z,y_1(:,1),'m')
legend('Is1')
xlabel('Lunghezza [m]')
ylabel('Densità di Potenza segnale 2[W/(m^2)]')
title('Andamento Intensità segnale 2')
grid on


figure
plot(z,y_2(:,1),'b')
legend('Is2')
xlabel('Lunghezza [m]')
ylabel('Densità di Potenza segnale 2[W/(m^2)]')
title('Andamento Intensità segnale 2')
grid on

figure
plot(z,y_0(:,2),'r')
legend('Ip0 con Ppo=1W')
xlabel('Lunghezza [m]')
ylabel('Densità di Potenza Pompa1  [W/(m^2)]')
title('Andamento Intensità Pompa 1')
grid on

figure
plot(z,y_1(:,2),'m')
legend('Ip1 con Pp1=2W')
xlabel('Lunghezza [m]')
ylabel('Densità di Potenza Pompa2 [W/(m^2)]')
title('Andamento Intensità Pompa 2')
grid on

figure
plot(z,y_2(:,2),'b')
legend('Ip2 con Pp2=3W')
xlabel('Lunghezza [m]')
ylabel('Densità di Potenza Pompa3 [W/(m^2)]')
title('Andamento Intensità Pompa 3')
grid on

G0=10*log10(y_0(:,1)/Is_0);
G1=10*log10(y_1(:,1)/Is_1);
G2=10*log10(y_2(:,1)/Is_2);
%Guadagno esatto derivante dalle slide
G_esatta1=10*log10(exp(g_r*Ip_0*z));

G_esatta2=10*log10(exp(g_r*Ip_1*z));

G_esatta3=10*log10(exp(g_r*Ip_2*z));

figure
plot(z,G0,'r',z,G_esatta1,'--')
legend('Guadagno ode45 P=1W','Guadagno esatto')
xlabel('Lunghezza [m]')
ylabel('Guadagno [dB]')
ylim([0 30])
title('Andamento Guadagno Raman')
grid on

figure
plot(z,G1,'m',z,G_esatta2,'--')
legend('Guadagno ode45 P=2W','Guadagno esatto')
xlabel('Lunghezza [m]')
ylabel('Guadagno [dB]')
ylim([0 30])
title('Andamento Guadagno Raman')
grid on

figure
plot(z,G2,'b',z,G_esatta3,'--')
legend('Guadagno ode45 P=3W','Guadagno esatto')
xlabel('Lunghezza [m]')
ylabel('Guadagno [dB]')
ylim([0 30])
title('Andamento Guadagno Raman')     
grid on


Codice: Seleziona tutto

function dy=Raman2(z,y)
%Co-propagating system equations
dy=zeros(2,1);

%y(1)==Is Intensity input Signal
%y(2)==Ip Intensity of the pump
%Is_dot=dy(1)=Is'=dIs/dz
%Ip_dot=dy(2)=Ip'=dIp/dz


g_r=10.^-13; %coeff. Raman [m/W]

Lambda_p=1450*10^(-9); %lunghezza d'onda pompa [nm]
Lambda_s=1550*10^(-9);%lunghezza d'onda segnale[nm]
fp=3e8/Lambda_p; % frequenza della pompa [Hz]
 fs=3e8/Lambda_s; %Frequenza del segnale [Hz]
 f_ratio=fp/fs;

alpha_s=0.046/1000; %attenuazione [dB/m]
alpha_p=alpha_s;
dy(1)=(g_r.*y(1).*y(2)-alpha_s.*y(1));
dy(2)=(-f_ratio.*g_r.*y(1).*y(2)-alpha_p.*y(2));

end

[dimaios]

Come no. A parte il segno di convoluzione anzichè la moltiplicazione e l'unità di misura sbagliata per l'ennesima volta.