Unable to reproduce a working code for threepointbvp

2 次查看（过去 30 天）
Dursman Mchabe
Dursman Mchabe 2018-8-28
评论： Dursman Mchabe 2018-8-28
Hi everyone,
I am struggling to reproduce a working threepointbvp code. Please see attached or below. What am I doing wrong?
The working code is:
function threebvp(solver)
%THREEBVP Three-point boundary value problem
% This example of multipoint BVPs comes from the study of a physiological
% flow problem described in C. Lin and L. Segel, Mathematics Applied to
% Deterministic Problems in the Natural Sciences, SIAM, Philadelphia, PA,
% 1988. For x in [0, lambda], the equations for the problem are
%
% v' = (C - 1)/n
% C' = (vC - min(x,1))/eta
%
% for known parameters n, kappa, and lambda > 1. eta = lambda^2/(n*kappa^2).
% The term min(x,1) in the equation for C'(x) is not smooth at x = 1.
% Lin and Segel describe this BVP as two problems, one set on [0, 1]
% and the other on [1, lambda], connected by the requirement that v(x)
% and C(x) be continuous at x = 1. The solution is to satisfy the boundary
% conditions v(0) = 0 and C(lambda) = 1. BVP4C solves this problem as
% a three-point BVP, imposing internal conditions at the interface point
% x = 1.
%
% This example solves the problem for n = 5e-2, lambda = 2, and a range
% kappa = 2,3,4,5. The solution for one value of kappa is used as guess
% for the next.
%
% By default, this example uses the BVP4C solver. Use syntax
% THREEBVP('bvp5c') to solve this problem with the BVP5C solver.
%
% See also BVP4C, BVP5C, BVPSET, BVPGET, BVPINIT, DEVAL, FUNCTION_HANDLE.
% Jacek Kierzenka and Lawrence F. Shampine
if nargin < 1
solver = 'bvp4c';
end
bvpsolver = fcnchk(solver);
% Problem parameter, shared with nested functions
n = 5e-2;
% Initial mesh - duplicate the interface point x = 1.
lambda = 2;
xinit = [0, 0.25, 0.5, 0.75, 1, 1, 1.25, 1.5, 1.75, 2];
% Constant initial guess for the solution
yinit = [1; 1];
% The initial profile
sol = bvpinit(xinit,yinit);
% For each kappa, the quantity of interest is the emergent osmolarity.
% This example compares the value computed by BVP4C, Os = 1/v(lambda),
% with Lin and Segel's asymptotic approximation.
fprintf(' kappa computed Os approximate Os \n')
for kappa = 2:5
eta = lambda^2/(n*kappa^2); % use in nested functions
sol = bvpsolver(@f,@bc,sol);
K2 = lambda*sinh(kappa/lambda)/(kappa*cosh(kappa));
approx = 1/(1 - K2);
computed = 1/sol.y(1,end);
fprintf(' %2i %10.3f %10.3f \n',kappa,computed,approx);
end
figure
plot(sol.x,sol.y(1,:),sol.x,sol.y(2,:),'--')
legend('v(x)','C(x)')
title(['A three-point BVP solved with ',upper(func2str(bvpsolver))])
xlabel(['\lambda = ',num2str(lambda),', \kappa = ',num2str(kappa),'.'])
ylabel('v and C')
% -----------------------------------------------------------------------
% Nested functions -- n and eta are provided by the outer function.
%
function dydx = f(x,y,region)
% Derivative function -- share n and eta with the outer function.
dydx = zeros(2,1);
dydx(1) = (y(2) - 1)/n;
% The definition of C'(x) depends on the region.
switch region
case 1 % x in [0 1]
dydx(2) = (y(1)*y(2) - x)/eta;
case 2 % x in [1 lambda]
dydx(2) = (y(1)*y(2) - 1)/eta;
otherwise
error('MATLAB:threebvp:BadRegionIndex','Incorrect region index: %d',region);
end
end
% -----------------------------------------------------------------------
function res = bc(YL,YR)
% Boundary (and internal) conditions
res = [ YL(1,1) % v(0) = 0
YR(1,1) - YL(1,2) % continuity of v(x) at x = 1
YR(2,1) - YL(2,2) % continuity of C(x) at x = 1
YR(2,end) - 1 ]; % C(lambda) = 1
end
% -----------------------------------------------------------------------
end % threebvp
where as my code is:
function threepointbvp(solver)
global A B C D E F G H I J K L M N O P Q R S T U V
%PARAMETERS
A = sym(9.4e-4); % ks = solid-liquid mass transfer (m/s)
B = sym(3.69e2); % Ap = solid-liquid surface area (m2/m3)
C = sym(2.89e-9); % DSO2 = diffusivity of SO2 (m2/s)
D = sym(1.6e-9); % DCO32 = diffusivity of CO32- (m2/s)
E = sym(4.92e-5); % CBs = concentration of Ca2+ on limestone surface (mol/m3)
F = sym(1.6e-8); % DH = diffusivity of H+ (m2/s)
G = sym(6.24); % KSO2 = SO2 dissociation equilibrium (mol/m3)
H = sym(2.35e-9); % DHSO3 = diffusivity of HSO3 (m2/s)
I = sym(5.68e-5); % KHSO3 = HSO3 dissociation equilibrium constant (mol/m3)
J = sym(1.69e-9); % DSO32 = diffusivity of SO32 (m2/s)
K = sym(149); % HSO2 = Henry's constant (m3Pa/mol)
L = sym(47.28325); % PSO2 = partial pressure of SO2 (Pa)
M = sym(4.14e-6); % kga = gas-side mass transfer product (1/s)
N = sym(8.314); % R = universal gas constant (J/mol.K)
O = sym(323.15); % T = reactor tempeature (K)
P = sym(2.09); % CH = concentration of H+ at the gas-liquid interface (mol/m3)
Q = sym(9.333E-06); % CH = concentration of H+ at the in the bulk slurry (mol/m3)
R = sym(0.679); % S(IV) = total sulfur concentration in the bulk slurry (mol/m3)
S = sym(4.88e-6); % C(IV) = total sulfur concentration in the bulk slurry (mol/m3)
T = sym(1.7e-3); % KCO2 = CO2 dissociation equilibrium constant (mol/m3)
U = 6.55e-8; % KHCO3 = HCO3 dissociation equilibrium constant (mol/m3)
V = 2.08e-9; % DHCO3 = Diffusivity of HCO3 (m2/s)
% VARIABLES
% y(1) = CSO2
% y(2) = CHSO3
% y(3) = CSO32
% y(4) = CHCO32
% y(5) = CCO32
% EQUATIONS
% Region I
% dy(1)/dx = y(6);
% d^2y(1)/dx^2 = (2*A*B / C)*(1 + (C*y(1) / 2*D*E) + (F*G*(y(2)/y(1)) / D*E)*E;
% dy(2)/dx = y(7);
% d^2y(2)/dx^2 = -(2*A*B / C)*(1 + (C*y(1) / 2*D*E) + (F*G*(y(2)/y(1)) /D*E)*E;
%Region II
% dy(3)/dx = y(8)
% dy(4)/dx = y(4);
% d^2y(3)/dx^2 = -(A*B / D)*(1 + (F*I*(y(3) / y(2)))/(D*E))*E;
% dy(5)/dx = y(5)
if nargin < 1
solver = 'bvp4c';
end
bvpsolver = fcnchk(solver);
xinit = [0.05, 0.10, 0.15, 0.20,0.25,0.30,0.35,0.40,0.45,0.50,0.5,0.55,0.60,0.65,0.70,0.75,0.80,0.85,0.90,0.95,1.00,];
yinit = [4.03e-08, 3.39e-02, 2.60e-01, 5.45e-01, 4.82e-03, 0, 0, 0];
solinit = bvpinit(xinit,yinit);
sol = bvpsolver(@materialbalanceodes,@materialbalancebc,solinit);
x = linspace(0.05,1,0.05);
y = deval(sol,x);
figure
plot(x,y(1,:), x,y(2,:), x,y(3,:), x,y(4,:), x,y(5,:), x,y(6,:), x,y(7,:), x,y(8,:))
legend('SO_2','HSO_{3}^{-}','SO_{3}^{2-}','HCO_{3}^{-}','CO_{3}^{2-}')
for x = 0.05:0.5
region = 'case 1';
for x = 0.55:1
region = 'case 2';
end
end
% -----------------------------------------------------------------------
function dydx = materialbalanceodes(x,y, region)
dydx = zeros(8,1);
switch region
case 1 % x in [0 0.5]
dydx(1) = y(6);
dydx(2) = (2 * A * B / C) * (1 + ((C * y(1)) / (2 * D * E)) + ((F * G )* (y(2)/y(1)) / (D*E)))*E;
dydx(3) = y(7);
dydx(4) = -((2 * A * B) / (H)) * (1 + ((C * y(1)) / (2 * D * E)) + ((F * G) * (y(2) / y(1)) / (D * E))) * E;
case 2 % x in [0.5 1]
dydx(5) = y(8);
dydx(6) = y(4);
dydx(7) = -((A * B) / D) * (1 + (F * I * (y(3) / y(2)))/(D * E)) * E;
dydx(8) = y(5);
otherwise
error('MATLAB:threebvp:BadRegionIndex','Incorrect region index: %d',region);
end
% -----------------------------------------------------------------------------------------
function res = materialbalancebc(ya,yb)
res = [ (ya(1,0.05) - (K * L)) % x = 0.05
(ya(2,0.05) - ((I * K * L) / P))
(ya(6,0.05) + M * ((L /( N * O))- (K * L)))
ya(1,0.5) % x = 0.5
ya(2,0.5) - ((R*G*Q)/(Q.^2 + G * Q + G * I))
ya(3,0.5)
ya(4,0.5)
- C * ya(6,0.5) - V * ya(4,0.5) - J * yb(8,0.5)
- H * yb(7,0.5) + H * yb(7,0.5) - V * ya(4,0.5) - 2 * J * ya(8,0.5)
(yb(2,1) - ((R * G * Q)/(Q.^2 + G * Q + G * I))) % x = 1
(yb(3,1) - ((R * G * I)/(Q.^2 + G * Q + G * I)))
yb(4,1) - ((S * T * Q)/(Q.^2 + T * Q + T * U))]
end
% ----------------------------------------------------------------------------------------------
end
end
Where am I going wrong?
Copyright 1984-2014 The MathWorks, Inc.
  3 个评论
显示 1更早的评论
Rik
Rik 2018-8-28
If it is working, why are you changing it? Also, why are you using global variables? They should be avoided if possible, or else you should have long names that will avoid any name collision (which would cause bugs that are very difficult to solve).
Dursman Mchabe
Dursman Mchabe 2018-8-28
Hi Rik, thanks a lot for the response. The reason I'm changing it is that I must use my equations and boundary conditions.

请先登录，再进行评论。

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

回答（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