Solving a boundary value problem using bvp4c: forcing a particular form of solution

10 visualizaciones (últimos 30 días)
Hello,
I am looking to replicate some results from a paper, linked here.
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.
  1 comentario
Jaya
Jaya el 24 de Abr. de 2021
Editada: Jaya el 24 de Abr. de 2021
Were you able to solve it later? Because I am also stuck with something similar.

Iniciar sesión para comentar.

Respuestas (0)

Community Treasure Hunt

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

Start Hunting!

Translated by