1 view (last 30 days)

Hello,

The results are a matlab plot that should generate as the figure in the paper, but I am having trouble using bvp5c, as discussed in the paper. The problem is that I cannot get the same form of the equation. Is there a way to force the solution to look more like that in the paper?

I provide some screenshots, and the relevant code I have been using here:

The code I have come up with so far is here:

%% Simulating the Raman Amplifier as done in:

% https://www.osapublishing.org/jlt/abstract.cfm?uri=jlt-33-18-3773

% Based on example of bvp4c

%% Defining the system of differential equations

for i = 2:1:20

solinit = bvpinit(linspace(0,7e3,i),@(x)init_sol(x));

options = bvpset('Stats','on','RelTol',1e-9);

sol = bvp5c(@ode,@bc,solinit,options);

% The solution at the mesh points

x = sol.x;

y = sol.y;

figure(i);

clf reset

plot(x,y','*')

title('Example problem for MUSN')

ylabel('bvp4c and MUSN (*) solutions')

xlabel('x')

shg

end

function dydx = ode(x,y)

%EX1ODE ODE function for Example 1 of the BVP tutorial.

% The components of y correspond to the original variables

% as y(1) = P_r, y(2) = P_p, y(3) = P_sbs.

g_r = 4.4e-14;

epsilon = 0.55;

g_b = 4e-13;

A_eff = 5.2e-11;

alpha_p = 2e-4;

alpha_r = 3e-4;

lambda_r = 1651e-9;

lambda_p = 1540e-9;

dydx = [ epsilon*g_r*y(2)*y(1)/A_eff - g_b*y(1)*y(3)/A_eff - alpha_r*y(1)

lambda_r/lambda_p * epsilon*g_r*y(2)*y(1)/A_eff + lambda_r/lambda_p * epsilon*g_r*y(2)*y(3)/A_eff + alpha_p*y(2)

-epsilon*g_r*y(2)*y(3)/A_eff - g_b*y(1)*y(3)/A_eff + alpha_r*y(3) ];

end

function res = bc(ya,yb)

%EX1BC Boundary conditions for Example 1 of the BVP tutorial.

% RES = EX1BC(YA,YB) returns a column vector RES of the

% residual in the boundary conditions resulting from the

% approximations YA and YB to the solution at the ends of

% the interval [a b]. The BVP is solved when RES = 0.

% The components of y correspond to the original variables

% as y(1) = P_r, y(2) = P_p, y(3) = P_sbs.

P_p = 4;

P_r = 6.5e-3;

P_sbs = 1e-6;

res = [ ya(1) - P_r

yb(2) - P_p

yb(3) - P_sbs];

end

function v = init_sol(x)

%EX1INIT guess for Example 1 of the BVP tutorial.

v = [6.5e-3

4

1e-10 ];

end

The desired output should look like this:

Any help or discussion would be greatly appreciated.

Opportunities for recent engineering grads.

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

Start Hunting!
## 0 Comments

Sign in to comment.