Using MATLAB (779505), страница 65
Текст из файла (страница 65)
Note that the boundary conditions are expressed in terms of theflux f rather than ∂u ⁄ ∂x . Also, of the two coefficients, only p can depend on u .MATLAB Partial Differential Equation SolverThis section describes:• The PDE solver, pdepe• PDE solver basic syntax• Additional PDE solver argumentsThe PDE SolverThe MATLAB PDE solver, pdepe, solves initial-boundary value problems forsystems of parabolic and elliptic PDEs in the one space variable x and time t .There must be at least one parabolic equation in the system.The pdepe solver converts the PDEs to ODEs using a second-order accuratespatial discretization based on a set of nodes specified by the user. Thediscretization method is described in [8].
The time integration is done withode15s. The pdepe solver exploits the capabilities of ode15s for solving thedifferential-algebraic equations that arise when Equation 15-3 contains ellipticequations, and for handling Jacobians with a specified sparsity pattern. ode15schanges both the time step and the formula dynamically.After discretization, elliptic equations give rise to algebraic equations. If theelements of the initial conditions vector that correspond to elliptic equations15-84Partial Differential Equationsare not “consistent” with the discretization, pdepe tries to adjust them beforebeginning the time integration. For this reason, the solution returned for theinitial time may have a discretization error comparable to that at any othertime.
If the mesh is sufficiently fine, pdepe can find consistent initial conditionsclose to the given ones. If pdepe displays a message that it has difficulty findingconsistent initial conditions, try refining the mesh. No adjustment is necessaryfor elements of the initial conditions vector that correspond to parabolicequations.PDE Solver Basic SyntaxThe basic syntax of the solver issol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan)Note Correspondences given are to terms used in “Introduction to PDEProblems” on page 15-83.The input arguments are:mSpecifies the symmetry of the problem. m can be 0 = slab,1 = cylindrical, or 2 = spherical. It corresponds to m in Equation 15-3.pdefunFunction that defines the components of the PDE.
It computes theterms c , f , and s in Equation 15-3, and has the form[c,f,s] = pdefun(x,t,u,dudx)where x and t are scalars, and u and dudx are vectors thatapproximate the solution u and its partial derivative with respectto x . c, f, and s are column vectors. c stores the diagonal elementsof the matrix c .icfunFunction that evaluates the initial conditions.
It has the formu = icfun(x)When called with an argument x, icfun evaluates and returns theinitial values of the solution components at x in the column vector u.15-8515Differential EquationsbcfunFunction that evaluates the terms p and q of the boundaryconditions. It has the form[pl,ql,pr,qr] = bcfun(xl,ul,xr,ur,t)where ul is the approximate solution at the left boundary xl = aand ur is the approximate solution at the right boundary xr = b .pl and ql are column vectors corresponding to p and the diagonal ofq evaluated at xl. Similarly, pr and qr correspond to xr.
Whenm > 0 and a = 0 , boundedness of the solution near x = 0 requiresthat the flux f vanish at a = 0 . pdepe imposes this boundarycondition automatically and it ignores values returned in pl and ql.xmeshVector [x0, x1, ..., xn] specifying the points at which a numericalsolution is requested for every value in tspan. x0 and xn correspondto a and b , respectively.Second-order approximation to the solution is made on the meshspecified in xmesh. Generally, it is best to use closely spaced meshpoints where the solution changes rapidly. pdepe does not select themesh in x automatically.
You must provide an appropriate fixedmesh in xmesh. The cost depends strongly on the length of xmesh.When m > 0 , it is not necessary to use a fine mesh near x = 0 toaccount for the coordinate singularity.The elements of xmesh must satisfy x0 < x1 < ... < xn. The length ofxmesh must be ≥ 3.tspanVector [t0, t1, ..., tf] specifying the points at which a solution isrequested for every value in xmesh. t0 and tf correspond to t 0 andt f , respectively.pdepe performs the time integration with an ODE solver that selectsboth the time step and formula dynamically.
The solutions at thepoints specified in tspan are obtained using the natural continuousextension of the integration formulas. The elements of tspan merelyspecify where you want answers and the cost depends weakly on thelength of tspan.The elements of tspan must satisfy t0 < t1 < ... < tf. The length oftspan must be ≥ 3.15-86Partial Differential EquationsThe output argument sol is a three-dimensional array, such that:• sol(:,:,k) approximates component k of the solution u .• sol(i,:,k) approximates component k of the solution at time tspan(i) andmesh points xmesh(:).• sol(i,j,k) approximates component k of the solution at time tspan(i) andthe mesh point xmesh(j).Additional PDE Solver ArgumentsFor more advanced applications, you can also specify as input arguments solveroptions and additional parameters that are passed to the PDE functions.optionsStructure of optional parameters that change the defaultintegration properties.
This is the seventh input argument.sol = pdepe(m,pdefun,icfun,bcfun,...xmesh,tspan,options)See “Improving PDE Solver Performance” on page 15-93 formore information.p1,p2...Parameters that the solver passes to pdefun, icfun, and bcfun.sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan,...options,p1,p2...)The solver passes any input parameters that follow the optionsargument to pdefun, icfun, and bcfun every time it calls them.Use options = [] as a placeholder if you set no options. In thepdefun argument list, parameters follow x, t, u, and dudx.f = pdefun(x,t,u,dudx,p1,p2,...)In the icfun argument list, parameters follow x.res = icfun(x,p1,p2,...)In the bcfun argument list, parameters follow xl, ul, xr, ur, andt.res = bcfun(xl,ul,xr,ur,t,p1,p2,...)See the pdex3 demo for an example.15-8715Differential EquationsRepresenting PDE ProblemsThis section describes:• The process for solving PDE problems using the MATLAB solver, pdepe• Evaluating the solution at specific pointsExample: A Single PDEThis example illustrates the straightforward formulation, solution, andplotting of the solution of a single PDE22 ∂u∂ uπ ------ = --------2∂t∂xThis equation holds on an interval 0 ≤ x ≤ 1 for times t ≥ 0 .
At t = 0 , thesolution satisfies the initial conditionu ( x, 0 ) = sin πxAt x = 0 and x = 1 , the solution satisfies the boundary conditionsu ( 0, t ) = 0πe–t∂u+ ------ ( 1, t ) = 0∂xNote The demo pdex1 contains the complete code for this example. The demouses subfunctions to place all functions it requires in a single M-file. To runthe demo type pdex1 at the command line. See “PDE Solver Basic Syntax” onpage 15-85 for more information.1 Rewrite the PDE. Write the PDE in the form∂u ∂u∂u∂u–m ∂ m ------ x f x, t, u, ------- + s x, t, u, -------c x, t, u, ------- ------- = x∂x ∂t∂x∂ x∂xThis is the form shown in Equation 15-3 and expected by pdepe.
See“Introduction to PDE Problems” on page 15-83 for more information. For thisexample, the resulting equation is15-88Partial Differential Equations2 ∂u0 ∂0 ∂uπ ------ = x ------ x ------- + 0∂t∂ x∂xwith parameter m = 0 and the terms∂u2c x, t, u, ------- = π∂x∂u∂uf x, t, u, ------- = ------∂x∂x∂us x, t, u, ------- = 0∂x2 Code the PDE in MATLAB.
Once you rewrite the PDE in the form shown above(Equation 15-3) and identify the terms, you can code the PDE in a function thatpdepe can use. The function must be of the form[c,f,s] = pdefun(x,t,u,dudx)where c, f, and s correspond to the c , f , and s terms. The code below computesc, f, and s for the example problem.function [c,f,s] = pdex1pde(x,t,u,DuDx)c = pi^2;f = DuDx;s = 0;3 Code the Initial Conditions Function. You must code the initial conditions in aMATLAB function of the formu = icfun(x)The code below represents the initial conditions in the MATLAB functionpdex1ic.function u0 = pdex1ic(x)u0 = sin(pi*x);4 Code the Boundary Conditions Function.
You must also code the boundaryconditions in a MATLAB function of the form[pl,ql,pr,qr] = bcfun(xl,ul,xr,ur,t)15-8915Differential EquationsThe boundary conditions, written in the same form as Equation 15-5, are∂uu(0, t) + 0 ⋅ ------(0, t) = 0∂xat x = 0andπe–t∂u+ 1 ⋅ ------(1, t) = 0∂xat x = 1The code below evaluates the components p ( x, t, u ) and q ( x, t ) of the boundaryconditions in the MATLAB function pdex1bc.function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t)pl = ul;ql = 0;pr = pi * exp(-t);qr = 1;In the function pdex1bc, pl and ql correspond to the left boundary conditions( x = 0 ), and pr and qr correspond to the right boundary condition ( x = 1 ).5 Select Mesh Points for the Solution. Before you use the MATLAB PDE solver, youneed to specify the mesh points ( t, x ) at which you want pdepe to evaluate thesolution.
Specify the points as vectors t and x.The vectors t and x play different roles in the solver (see “MATLAB PartialDifferential Equation Solver” on page 15-84). In particular, the cost and theaccuracy of the solution depend strongly on the length of the vector x. However,the computation is much less sensitive to the values in the vector t.This example requests the solution on the mesh produced by 20 equally spacedpoints from the spatial interval [0,1] and five values of t from the time interval[0,2].x = linspace(0,1,20);t = linspace(0,2,5);6 Apply the PDE Solver.
The example calls pdepe with m = 0, the functionspdex1pde, pdex1ic, and pdex1bc, and the mesh defined by x and t at whichpdepe is to evaluate the solution. The pdepe function returns the numericalsolution in a three-dimensional array sol, where sol(i,j,k) approximates thekth component of the solution, u k , evaluated at t(i) and x(j).15-90Partial Differential Equationsm = 0;sol = pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t);This example uses @ to pass pdex1pde, pdex1ic, and pdex1bc as functionhandles to pdepe.Note See the function_handle (@), func2str, and str2func reference pages,and the “Function Handles” chapter of “Programming and Data Types” in theMATLAB documentation for information about function handles.7 View the Results. Complete the example by displaying the results:• Extract and display the first solution component.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
