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

3 次查看（过去 30 天）

显示更早的评论

Michael Soskind 2020-5-7

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/523778-solving-a-boundary-value-problem-using-bvp4c-forcing-a-particular-form-of-solution

编辑： Jaya 2021-4-24

在 MATLAB Online 中打开

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