role of flux f in bcfun of pdepe for systems of PDEs
3 次查看(过去 30 天)
显示 更早的评论
feynman feynman
2024-3-3
To solve utt=uxx, as the flux term f has a 0 entry in pdefun, the same f is in bcfun. In bcfun, there's q(x,t)f(x,t,u,ux), which appears to mean that when f=[0;ux(1)], q=[1;0] is identical to q=[0;0] because the 1st entry of q doesn't matter after multiplying the 0 in f. But when q=[1;0] is used, the program yields wrong results, why?
wave()
function wave
x=linspace(-pi,pi);t=linspace(0,2*pi);
sol=pdepe(0,@pdefun,@icfun,@bcfun,x,t);
u=sol(:,:,1);
surf(x,t,u)
function [c,f,s]=pdefun(x,t,u,ux)
c=[1;1];
f=[0;ux(1)];
s=[u(2);0];
end
function u0=icfun(x)
u0=[sin(x);0];
end
function [pl,ql,pr,qr]=bcfun(xl,ul,xr,ur,t)
% p(x,t,u)+q(x,t)f(x,t,u,ux)=0
pl=ul;pr=ur;ql=[1;0];qr=[1;0];
end
end
2 个评论
Torsten
2024-3-3
If you try to solve the standard-hyperbolic equation (wave equation) with a solver for parabolic-elliptic equations (for which f should always be different from 0), you shouldn't be surprised to get results that cannot be explained.
采纳的回答
Torsten
2024-3-4
移动:Torsten
2024-3-6
This multiplication is nowhere performed within "pdepe". But if you look at the pdepe code below, the cases where q = 0 and q different from 0 are treated differently (independent of f).
type pdepe
function varargout = pdepe(m,pde,ic,bc,xmesh,t,options,varargin)
%PDEPE Solve initial-boundary value problems for parabolic-elliptic PDEs in 1-D.
% SOL = PDEPE(M,PDEFUN,ICFUN,BCFUN,XMESH,TSPAN) solves initial-boundary
% value problems for small systems of parabolic and elliptic PDEs in one
% space variable x and time t to modest accuracy. There are npde unknown
% solution components that satisfy a system of npde equations of the form
%
% c(x,t,u,Du/Dx) * Du/Dt = x^(-m) * D(x^m * f(x,t,u,Du/Dx))/Dx + s(x,t,u,Du/Dx)
%
% Here f(x,t,u,Du/Dx) is a flux and s(x,t,u,Du/Dx) is a source term. m must
% be 0, 1, or 2, corresponding to slab, cylindrical, or spherical symmetry,
% respectively. The coupling of the partial derivatives with respect to
% time is restricted to multiplication by a diagonal matrix c(x,t,u,Du/Dx).
% The diagonal elements of c are either identically zero or positive.
% An entry that is identically zero corresponds to an elliptic equation and
% otherwise to a parabolic equation. There must be at least one parabolic
% equation. An entry of c corresponding to a parabolic equation is permitted
% to vanish at isolated values of x provided they are included in the mesh
% XMESH, and in particular, is always allowed to vanish at the ends of the
% interval. The PDEs hold for t0 <= t <= tf and a <= x <= b. The interval
% [a,b] must be finite. If m > 0, it is required that 0 <= a. The solution
% components are to have known values at the initial time t = t0, the
% initial conditions. The solution components are to satisfy boundary
% conditions at x=a and x=b for all t of the form
%
% p(x,t,u) + q(x,t) * f(x,t,u,Du/Dx) = 0
%
% q(x,t) is a diagonal matrix. The diagonal elements of q must be either
% identically zero or never zero. Note that the boundary conditions are
% expressed in terms of the flux rather than Du/Dx. Also, of the two
% coefficients, only p can depend on u.
%
% The input argument M defines the symmetry of the problem. PDEFUN, ICFUN,
% and BCFUN are function handles.
%
% [C,F,S] = PDEFUN(X,T,U,DUDX) evaluates the quantities defining the
% differential equation. The input arguments are scalars X and T and
% vectors U and DUDX that approximate the solution and its partial
% derivative with respect to x, respectively. PDEFUN returns column
% vectors: C (containing the diagonal of the matrix c(x,t,u,Dx/Du)),
% F, and S (representing the flux and source term, respectively).
%
% U = ICFUN(X) evaluates the initial conditions. For a scalar X, ICFUN
% must return a column vector, corresponding to the initial values of
% the solution components at X.
%
% [PL,QL,PR,QR] = BCFUN(XL,UL,XR,UR,T) evaluates the components of the
% boundary conditions at time T. XL and XR are scalars representing the
% left and right boundary points. UL and UR are column vectors with the
% solution at the left and right boundary, respectively. PL and QL are
% column vectors corresponding to p and the diagonal of q, evaluated at
% the left boundary, similarly PR and QR correspond to the right boundary.
% When m > 0 and a = 0, boundedness of the solution near x = 0 requires
% that the flux f vanish at a = 0. PDEPE imposes this boundary condition
% automatically.
%
% PDEPE returns values of the solution on a mesh provided as the input
% array XMESH. The entries of XMESH must satisfy
% a = XMESH(1) < XMESH(2) < ... < XMESH(NX) = b
% for some NX >= 3. Discontinuities in c and/or s due to material
% interfaces are permitted if the problem requires the flux f to be
% continuous at the interfaces and a mesh point is placed at each
% interface. The ODEs resulting from discretization in space are integrated
% to obtain approximate solutions at times specified in the input array
% TSPAN. The entries of TSPAN must satisfy
% t0 = TSPAN(1) < TSPAN(2) < ... < TSPAN(NT) = tf
% for some NT >= 3. The arrays XMESH and TSPAN do not play the same roles
% in PDEPE: The time integration is done with an ODE solver that selects
% both the time step and formula dynamically. The cost depends weakly on
% the length of TSPAN. Second order approximations to the solution are made
% on the mesh specified in XMESH. Generally it is best to use closely
% spaced points where the solution changes rapidly. PDEPE does not select
% the mesh in x automatically like it does in t; you must choose an
% appropriate fixed mesh yourself. The discretization takes into account
% the coordinate singularity at x = 0 when m > 0, so it is not necessary to
% use a fine mesh near x = 0 for this reason. The cost depends strongly on
% the length of XMESH.
%
% The solution is returned as a multidimensional array SOL. UI = SOL(:,:,i)
% is an approximation to component i of the solution vector u for
% i = 1:npde. The entry UI(j,k) = SOL(j,k,i) approximates UI at
% (t,x) = (TSPAN(j),XMESH(k)).
%
% SOL = PDEPE(M,PDEFUN,ICFUN,BCFUN,XMESH,TSPAN,OPTIONS) solves as above
% with default integration parameters replaced by values in OPTIONS, an
% argument created with the ODESET function. Only some of the options of
% the underlying ODE solver are available in PDEPE - RelTol, AbsTol,
% NormControl, InitialStep, and MaxStep. See ODESET for details.
%
% [SOL,TSOL,SOLE,TE,IE] = PDEPE(M,PDEFUN,ICFUN,BCFUN,XMESH,TSPAN,OPTIONS)
% with the 'Events' property in OPTIONS set to a function handle EVENTS,
% solves as above while also finding where event functions g(t,u(x,t))
% are zero. For each function you specify whether the integration is to
% terminate at a zero and whether the direction of the zero crossing
% matters. These are the three column vectors returned by EVENTS:
% [VALUE,ISTERMINAL,DIRECTION] = EVENTS(M,T,XMESH,UMESH).
% XMESH contains the spatial mesh and UMESH is the solution at the mesh
% points. Use PDEVAL to evaluate the solution between mesh points.
% For the I-th event function: VALUE(I) is the value of the function,
% ISTERMINAL(I) = 1 if the integration is to terminate at a zero of this
% event function and 0 otherwise. DIRECTION(I) = 0 if all zeros are to be
% computed (the default), +1 if only zeros where the event function is
% increasing, and -1 if only zeros where the event function is decreasing.
% Output TSOL is a column vector of times specified in TSPAN, prior to
% first terminal event. SOL(j,:,:) is the solution at T(j). TE is a vector
% of times at which events occur. SOLE(j,:,:) is the solution at TE(j) and
% indices in vector IE specify which event occurred.
%
% If UI = SOL(j,:,i) approximates component i of the solution at time TSPAN(j)
% and mesh points XMESH, PDEVAL evaluates the approximation and its partial
% derivative Dui/Dx at the array of points XOUT and returns them in UOUT
% and DUOUTDX: [UOUT,DUOUTDX] = PDEVAL(M,XMESH,UI,XOUT)
% NOTE that the partial derivative Dui/Dx is evaluated here rather than the
% flux. The flux is continuous, but at a material interface the partial
% derivative may have a jump.
%
% Example
% x = linspace(0,1,20);
% t = [0 0.5 1 1.5 2];
% sol = pdepe(0,@pdex1pde,@pdex1ic,@pdex1bc,x,t);
% solve the problem defined by the function pdex1pde with the initial and
% boundary conditions provided by the functions pdex1ic and pdex1bc,
% respectively. The solution is obtained on a spatial mesh of 20 equally
% spaced points in [0 1] and it is returned at times t = [0 0.5 1 1.5 2].
% Often a good way to study a solution is plot it as a surface and use
% Rotate 3D. The first unknown, u1, is extracted from sol and plotted with
% u1 = sol(:,:,1);
% surf(x,t,u1);
% PDEX1 shows how this problem can be coded using subfunctions. For more
% examples see PDEX2, PDEX3, PDEX4, PDEX5. The examples can be read
% separately, but read in order they form a mini-tutorial on using PDEPE.
%
% See also PDEVAL, ODE15S, ODESET, ODEGET, FUNCTION_HANDLE.
% The spatial discretization used is that of R.D. Skeel and M. Berzins, A
% method for the spatial discretization of parabolic equations in one space
% variable, SIAM J. Sci. Stat. Comput., 11 (1990) 1-32. The time
% integration is done with ODE15S. PDEPE exploits the capabilities this
% code has for solving the differential-algebraic equations that arise when
% there are elliptic equations and for handling Jacobians with specified
% sparsity pattern.
% Elliptic equations give rise to algebraic equations after discretization.
% Initial values taken from the given initial conditions may not be
% "consistent" with the discretization, so PDEPE may have to adjust these
% values slightly before beginning the time integration. For this reason
% the values output for these components at t0 may have a discretization
% error comparable to that at any other time. No adjustment is necessary
% for solution components corresponding to parabolic equations. If the mesh
% is sufficiently fine, the program will be able to find consistent initial
% values close to the given ones, so if it has difficulty initializing,
% try refining the mesh.
% Lawrence F. Shampine and Jacek Kierzenka
% Copyright 1984-2022 The MathWorks, Inc.
% $Revision: 1.8.4.9.12.1 $ $Date: 2013/09/27 03:10:20 $
if nargin < 7
options = [];
end
switch m
case {0, 1, 2}
otherwise
error(message('MATLAB:pdepe:InvalidM'))
end
nt = length(t);
if nt < 3
error(message('MATLAB:pdepe:TSPANnotEnoughPts'))
end
if any(diff(t) <= 0)
error(message('MATLAB:pdepe:TSPANnotIncreasing'))
end
xmesh = xmesh(:);
if m > 0 && xmesh(1) < 0
error(message('MATLAB:pdepe:NegXMESHwithPosM'))
end
nx = length(xmesh);
if nx < 3
error(message('MATLAB:pdepe:XMESHnotEnoughPts'))
end
if any(diff(xmesh) <= 0)
error(message('MATLAB:pdepe:XMESHnotIncreasing'))
end
% Initialize the nx-1 points xi where functions will be evaluated
% and as many coefficients in the difference formulas as possible.
% For problems with coordinate singularity, use 'singular' formula
% on all subintervals.
singular = (xmesh(1) == 0 && m > 0);
xL = xmesh(1:end-1);
xR = xmesh(2:end);
xM = xL + 0.5*(xR-xL);
switch m
case 0
xi = xM;
zeta = xi;
case 1
if singular
xi = (2/3)*(xL.^2 + xL.*xR + xR.^2) ./ (xL+xR);
else
xi = (xR-xL) ./ log(xR./xL);
end
zeta = (xi .* xM).^(1/2);
case 2
if singular
xi = (2/3)*(xL.^2 + xL.*xR + xR.^2) ./ (xL+xR);
else
xi = xL .* xR .* log(xR./xL) ./ (xR-xL);
end
zeta = (xL .* xR .* xM).^(1/3);
end
if singular
xim = (zeta .^ (m+1))./xi;
else
xim = xi .^ m;
end
zxmp1 = (zeta.^(m+1) - xL.^(m+1)) / (m+1);
xzmp1 = (xR.^(m+1) - zeta.^(m+1)) / (m+1);
% Form the initial values with a column of unknowns at each
% mesh point and then reshape into a column vector for dae.
temp = feval(ic,xmesh(1),varargin{:});
if size(temp,1) ~= numel(temp)
error(message('MATLAB:pdepe:InvalidOutputICFUN'))
end
npde = length(temp);
y0 = zeros(npde,nx);
y0(:,1) = temp;
for j = 2:nx
y0(:,j) = feval(ic,xmesh(j),varargin{:});
end
% Classify the equations so that a constant, diagonal mass matrix
% can be formed. The entries of c are to be identically zero or
% vanish only at the entries of xmesh, so need be checked only at one
% point not in the mesh.
[U,Ux] = pdentrp(singular,m,xmesh(1),y0(:,1),xmesh(2),y0(:,2),xi(1));
[c,f,s] = feval(pde,xi(1),t(1),U,Ux,varargin{:});
if any([size(c,1),size(f,1),size(s,1)]~=npde) || ...
any([numel(c),numel(f),numel(s)]~=npde)
error(message('MATLAB:pdepe:UnexpectedOutputPDEFUN',sprintf('%d',npde)))
end
[pL,qL,pR,qR] = feval(bc,xmesh(1),y0(:,1),xmesh(nx),y0(:,nx),t(1),varargin{:});
if any([size(pL,1),size(qL,1),size(pR,1),size(qR,1)]~=npde) || ...
any([numel(pL),numel(qL),numel(pR),numel(qR)]~=npde)
error(message('MATLAB:pdepe:UnexpectedOutputBCFUN',sprintf('%d',npde)))
end
D = ones(npde,nx);
D( c == 0, 2:nx-1) = 0;
if ~singular
D( qL == 0, 1) = 0;
end
D( qR == 0, nx) = 0;
M = spdiags( D(:), 0, npde*nx, npde*nx);
% Construct block-diagonal pattern of Jacobian matrix
S = kron( spdiags(ones(nx,3),-1:1,nx,nx), ones(npde,npde));
% Extract relevant options and augment with new ones
reltol = odeget(options,'RelTol',[]);
abstol = odeget(options,'AbsTol',[]);
normcontrol = odeget(options,'NormControl',[]);
initialstep = odeget(options,'InitialStep',[]);
maxstep = odeget(options,'MaxStep',[]);
events = odeget(options,'Events',[]); % events(m,t,xmesh,umesh)
hasEvents = ~isempty(events);
if hasEvents
eventfcn = @(t,y) events(m,t,xmesh,y,varargin{:});
else
eventfcn = [];
end
opts = odeset('RelTol',reltol,'AbsTol',abstol,'NormControl',normcontrol,...
'InitialStep',initialstep,'MaxStep',maxstep,'Events',eventfcn,...
'Mass',M,'JPattern',S);
% Call DAE solver
tfinal = t(end);
try
if hasEvents
[t,y,te,ye,ie] = ode15s(@pdeodes,t,y0(:),opts);
else
[t,y] = ode15s(@pdeodes,t,y0(:),opts);
end
catch ME
if strcmp(ME.identifier,'MATLAB:daeic12:IndexGTOne')
error(message('MATLAB:pdepe:SpatialDiscretizationFailed'))
else
rethrow(ME);
end
end
% Verify and process the solution
if t(end) ~= tfinal
nt = length(t);
if ~hasEvents || (te(end) ~= t(end)) % did not stop on a terminal event
warning(message('MATLAB:pdepe:TimeIntegrationFailed',sprintf('%e',t(end))));
end
end
sol = zeros(nt,nx,npde);
for i = 1:npde
sol(:,:,i) = y(:,i:npde:end);
end
varargout{1} = sol;
if hasEvents % [sol,t,sole,te,ie] = pdepe(...)
varargout{2} = t(:);
sole = zeros(length(te),nx,npde);
for i = 1:npde
sole(:,:,i) = ye(:,i:npde:end);
end
varargout{3} = sole;
varargout{4} = te(:);
varargout{5} = ie(:);
end
%---------------------------------------------------------------------------
% Nested functions
%---------------------------------------------------------------------------
function dudt = pdeodes(tnow,y)
%PDEODES Assemble the difference equations and evaluate the time derivative
% for the ODE system.
u = reshape(y,npde,nx);
up = zeros(npde,nx);
[U,Ux] = pdentrp(singular,m,xmesh(1),u(:,1),xmesh(2),u(:,2),xi(1));
[cL,fL,sL] = feval(pde,xi(1),tnow,U,Ux,varargin{:});
% Evaluate the boundary conditions
[pL,qL,pR,qR] = feval(bc,xmesh(1),u(:,1),xmesh(nx),u(:,nx),tnow,varargin{:});
% Left boundary
if singular
denom = cL;
denom(denom == 0) = 1;
up(:,1) = (sL + (m+1) * fL / xi(1)) ./ denom;
else
up(:,1) = pL;
idx = (qL ~= 0);
denom = (qL(idx)/xmesh(1)^m) .* (zxmp1(1)*cL(idx));
denom(denom == 0) = 1;
up(idx,1) = ( pL(idx) + (qL(idx)/xmesh(1)^m) .* ...
(xim(1)*fL(idx) + zxmp1(1)*sL(idx))) ./ denom;
end
% Interior points
for ii = 2:nx-1
[U,Ux] = pdentrp(singular,m,xmesh(ii),u(:,ii),xmesh(ii+1),u(:,ii+1),xi(ii));
[cR,fR,sR] = feval(pde,xi(ii),tnow,U,Ux,varargin{:});
denom = zxmp1(ii) * cR + xzmp1(ii-1) * cL;
denom(denom == 0) = 1;
up(:,ii) = ((xim(ii) * fR - xim(ii-1) * fL) + ...
(zxmp1(ii) * sR + xzmp1(ii-1) * sL)) ./ denom;
cL = cR;
fL = fR;
sL = sR;
end
% Right boundary
up(:,nx) = pR;
idx = (qR ~= 0);
denom = -(qR(idx)/xmesh(nx)^m) .* (xzmp1(nx-1)*cL(idx));
denom(denom == 0) = 1;
up(idx,nx) = ( pR(idx) + (qR(idx)/xmesh(nx)^m) .* ...
(xim(nx-1)*fL(idx) - xzmp1(nx-1)*sL(idx))) ./ denom;
dudt = up(:);
end % pdeodes
% --------------------------------------------------------------------------
end % pdepe
7 个评论
feynman feynman
2024-3-4
移动:Torsten
2024-3-6
Oh! Then why is such a misleading q(x,t)f(x,t,u,ux) put in the documentation? Also do q=1 or q=2 matter?
Torsten
2024-3-4
移动:Torsten
2024-3-6
Analyze this part of the code and your questions will be answered:
% Evaluate the boundary conditions
[pL,qL,pR,qR] = feval(bc,xmesh(1),u(:,1),xmesh(nx),u(:,nx),tnow,varargin{:});
% Left boundary
if singular
denom = cL;
denom(denom == 0) = 1;
up(:,1) = (sL + (m+1) * fL / xi(1)) ./ denom;
else
up(:,1) = pL;
idx = (qL ~= 0);
denom = (qL(idx)/xmesh(1)^m) .* (zxmp1(1)*cL(idx));
denom(denom == 0) = 1;
up(idx,1) = ( pL(idx) + (qL(idx)/xmesh(1)^m) .* ...
(xim(1)*fL(idx) + zxmp1(1)*sL(idx))) ./ denom;
end
% Interior points
for ii = 2:nx-1
[U,Ux] = pdentrp(singular,m,xmesh(ii),u(:,ii),xmesh(ii+1),u(:,ii+1),xi(ii));
[cR,fR,sR] = feval(pde,xi(ii),tnow,U,Ux,varargin{:});
denom = zxmp1(ii) * cR + xzmp1(ii-1) * cL;
denom(denom == 0) = 1;
up(:,ii) = ((xim(ii) * fR - xim(ii-1) * fL) + ...
(zxmp1(ii) * sR + xzmp1(ii-1) * sL)) ./ denom;
cL = cR;
fL = fR;
sL = sR;
end
% Right boundary
up(:,nx) = pR;
idx = (qR ~= 0);
denom = -(qR(idx)/xmesh(nx)^m) .* (xzmp1(nx-1)*cL(idx));
denom(denom == 0) = 1;
up(idx,nx) = ( pR(idx) + (qR(idx)/xmesh(nx)^m) .* ...
(xim(nx-1)*fL(idx) - xzmp1(nx-1)*sL(idx))) ./ denom;
feynman feynman
2024-3-5
移动:Torsten
2024-3-6
Thank you. Though this tells a lot I don't know what this is doing:
up(idx,1) = ( pL(idx) + (qL(idx)/xmesh(1)^m) .* ...
(xim(1)*fL(idx) + zxmp1(1)*sL(idx))) ./ denom;
and how it relates to q(x,t)f(x,t,u,ux).
Torsten
2024-3-5
移动:Torsten
2024-3-6
If q = 0, the values of u in the boundary points are set to p:
up(:,1) = pL;
up(:,nx) = pR;
If q is not equal to 0 in equation idx, the values of u(idx) in the boundary points are set to
up(idx,1) = ( pL(idx) + (qL(idx)/xmesh(1)^m) .* ...
(xim(1)*fL(idx) + zxmp1(1)*sL(idx))) ./ denom;
up(idx,nx) = ( pR(idx) + (qR(idx)/xmesh(nx)^m) .* ...
(xim(nx-1)*fL(idx) - xzmp1(nx-1)*sL(idx))) ./ denom;
Thus the product f*q is nowhere checked or used.
And if you now ask next about why p+q*f = 0 leads to the above updates for u in the boundary points: I don't know.
But it should all be available in the Berzins publication.
feynman feynman
2024-3-5
移动:Torsten
2024-3-6
Your penultimate sentence is exactly what I don't understand even after reading Berzins paper :)
Torsten
2024-3-5
移动:Torsten
2024-3-6
But your question is answered.
As far as I remember, you asked why pdepe handles the cases q = [0 0] and q = [1 0] differently. And this happens because up(1) is computed differently.
If you have a new question, you ahould accept the answer given here and open a new post.
By the way: Paragraph 2.3 of the paper explicitly explains how the equations in the boundary points are derived.
feynman feynman
2024-3-6
移动:Torsten
2024-3-6
I appreciate your help a lot, but you put all answers in comments so I can't accept your answer by clicking any button.
更多回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Eigenvalue Problems 的更多信息
标签
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!发生错误
由于页面发生更改,无法完成操作。请重新加载页面以查看其更新后的状态。
您也可以从以下列表中选择网站:
如何获得最佳网站性能
选择中国网站(中文或英文)以获得最佳网站性能。其他 MathWorks 国家/地区网站并未针对您所在位置的访问进行优化。
美洲
- América Latina (Español)
- Canada (English)
- United States (English)
欧洲
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
亚太
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)