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

1 次查看（过去 30 天）
Michael Soskind
Michael Soskind 2020-5-7
编辑： Jaya 2021-4-24
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.

请先登录，再回答此问题。

回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Boundary Value Problems 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by