ODEs system

i'm trying to solve this system of differential equations, can someone told me how to use the ODE45 function?

dPp/dz= Pp*gammap*(sigmape*N2 −sigmapa*N1)−α*Pp;

dPs/dz= Ps*gammas*(sigmase*N2 −sigmasa*N1)−α*Ps:

dPASE/dz=PASE*gammas*(sigmaSE*N2-sigmasa*N1)+2*sigmase*h*gammas*Vs*Δv-alfas*PASE;

with N1=ρ*(1+W12*t)/(1+(W12+W21)*t+R*t)

N2=ρ*(R*t+W21*t)/(1+(W12+W21)*t+R*t)

W12=[(sigmasa*gammas) / (h*Vs*A)](Ps+PASE) W21=[(sigmase*gammas) / (h*Vs*A)](Ps+PASE) R=[(Pp*gammap*sigmapa) / (h*Vp*A)](Ps+PASE)

(gammap,gammas,sigmase,sigmape,sigmapa,sigmasa,h,Vs,Vp,A,Δv,α,ρ are known parameters).

Answer 1

[Arnaud Miege]

Also, there seems to be 4 variables, not three. Below is the modified code I used and the results I get with your parameters. Again, check the values and the units.

Pp=y(1);
Ps=y(2);
PASE_plus=y(3);
PASE_minus=y(4);

W12=(((sigmasa*gammas)/(h*Vs*A))*(Ps+PASE_plus+PASE_minus));
W21=(((sigmase*gammas)/(h*Vs*A))*(Ps+PASE_plus+PASE_minus));
R=((Pp*gammap*sigmapa)/(h*Vp*A));
N1=(rho*((1+W21*t)/(1+(W12+W21)*t+R*t)));
N2=(rho*((R*t+W12*t)/(1+(W12+W21)*t+R*t)));

dydz(1)=Pp*gammap*(sigmape*N2-sigmapa*N1)-alfap*Pp;
dydz(2)=Ps*gammas*(sigmase*N2-sigmasa*N1)-alfas*Ps;
dydz(3)=PASE_plus*gammas*(sigmase*N2- ... sigmasa*N1)+2*sigmase*N2*gammas*h*Vs*deltav-alfas*PASE_plus;
dydz(4)=-PASE_minus*gammas*(sigmase*N2- ... sigmasa*N1)+2*sigmase*N2*gammas*h*Vs*deltav+alfas*PASE_minus;

And this is how I call the|ode| solver:

[z,y]=ode45('signalFW',[0 20],[10 0.001 0 0]);
figure(1),subplot(2,1,1),plot(z,y(:,1:2)),grid on,ylabel('Power in mW'),legend('Pp','Ps');
subplot(2,1,2),plot(z,y(:,3:4)),grid on,xlabel('length EDFA in m'),ylabel('Power in mW'),legend('PASE+','PASE-');

This gives me the following results:

Is that what you'd expect? If not, check the values of your parameters.

Arnaud

[Andrea]

Your results are the same as mine..maybe there is something wrong with the parameters. Thank you!

Answer 2

[Arnaud Miege]

Have a look at the examples provided in the documentation. You need to write a function that computes the derivatives of your variables as a function of time and the variables themselves, and pass this to the ode solver.

HTH,

Arnaud

[Andrea]

Thanks you for the answer.

I've built a m-file for the derivatives:

function dydz=signalFW(z,y)

dydz=zeros(size(y));

sigmapa=2.91e-25;

sigmasa=2.86e-25;

sigmase=4.03e-25;

sigmape=0.32e-25;

rho=4.4e24;

gammas=0.2;

gammap=0.25;

h=6.63e-34;

Vs=(3e8)/(1550e-9);

Vp=(3e8)/(1480e-9);

A=1.98e-11;

t=0.011;

alfas=4.6e-5;

alfap=5.52e-5;

deltav=3100e9;

Pp=y(1);

Ps=y(2);

PASEF=y(3);

W12=(((sigmasa*gammas)/(h*Vs*A))*(Ps+PASEF));

W21=(((sigmase*gammas)/(h*Vs*A))*(Ps+PASEF));

R=((Pp*gammap*sigmapa)/(h*Vp*A));

N1=(rho*((1+W21*t)/(1+(W12+W21)*t+R*t)));

N2=(rho*((R*t+W12*t)/(1+(W12+W21)*t+R*t)));

dydz(1)=Pp*gammap*(sigmape*N2-sigmapa*N1)-alfap*Pp;

dydz(2)=Ps*gammas*(sigmase*N2-sigmasa*N1)-alfas*Ps;

dydz(3)=PASEF*gammas*(sigmase*N2-

sigmasa*N1)+2*sigmase*N2*gammas*h*Vs*deltav-alfas*PASEF;

and I've passed it in this ODE file:

[z,y]=ode45('signalFWprova',[0 20],[10 0.001 0]);

figure(1),plot(z,y),xlabel('length EDFA in m'),ylabel('Power in mW'),legend('Pp','Ps','PASEF');

The result are not so good. Are this programs correct?

[Arnaud Miege]

I noticed there's no dependency on z in the calculation of dydz. Is that intentional? I can't really comment on the correctness of your equations, it's up to you to check that. The solver is able to solve the equations, so if the results are not what you expect, you probably want to double-check your equations are correct.

I also notice that you are multiplying/dividing some pretty small numbers by some pretty large ones. You might want to try and scale some of your parameters/variables to get a better conditioned system.

[Andrea]

Ps, Pp and PASE are signals which evolve in space. Do I have to write Ps(z), Pp(z) and PASE(z) in the differential equations? Otherwise what can I write?

P.S. the equations are correct, i've double-checked them

[Arnaud Miege]

If the derivatives of these variables do not depend on z through the differential equations (which is what you have in your code), then that's fine. If they do depend on z, then something's not right. The equations may be right, but have you checked they're coded correctly in MATLAB? Also, what do you mean by "the results are not so good".

[Andrea]

The result are not so good because in the original article (where i found the system), they are completely different. I'm not sure that my matlab code is correct...If the derivatives depend on z what i have to change in my program? (i think that error is probably in the code). If you want i can send you the article

[Matt Tearle]

Minor dumb thing: your function is signalFW but you're calling ode45 on signalFWprova. Assuming these are the same, I get sensible-looking results, so maybe you should post a link to the article or something like that so we can see what the results should look like.

[Arnaud Miege]

Good point, thanks Matt.

[Andrea]

The article is: "Modeling, optimization, and experimental evaluation of remotely pumped double-pass EDFA" and you can find it here http://onlinelibrary.wiley.com/doi/10.1002/mop.22687/abstract (you need an account to see it).

There is also: "A computer based simulator for Erbium-Doped Fiber Amplifier" here http://ieeexplore.ieee.org/xpls/abs_all.jsp?arnumber=5556859. The correct model is in the second article while the result are in the first

[Arnaud Miege]

I don't have an account, so can't access either publication.

[Andrea]

Can I send you an e-mail?

[Arnaud Miege]

Can you just upload an image of the equations and the results you expect for certain numerical parameters?

http://www.mathworks.com/matlabcentral/answers/7924-where-can-i-upload-images-and-files-for-use-on-matlab-answers

Answer 3

[Arnaud Miege]

This is the results I get with your code:

[Andrea]

You can find the model and the result in these two articles (the correct model is in the second)

http://www.sendspace.com/file/uujskz

http://www.sendspace.com/file/c7abb0

Is my code correct (for the first iteration of the article's model?)

[Arnaud Miege]

The code looks correct if you refer to equations (1) to (9) of the paper. However, you should check that the units used for the various parameters are consistent, and that the numerical values make sense.