Using MATLAB (779505), страница 63
Текст из файла (страница 63)
Torun this example type shockbvp at the command line. See “BVP Solver BasicSyntax” on page 15-61 and “Representing BVP Problems” on page 15-63 formore information.Note This problem appears in [1] to illustrate the mesh selection capabilityof a well established BVP code COLSYS.1 Code the ODE and Boundary Condition Functions. Code the differential equationand the boundary conditions as functions that bvp4c can use. Because there isan additional known parameter ε , the functions must be of the formdydx = odefun(x,y,p1)res = bcfun(ya,yb,p1)The code below represents the differential equation and the boundaryconditions in the MATLAB functions shockODE and shockBC. Note thatshockODE is vectorized to improve solver performance.
The additionalparameter ε is represented by e.function dydx = shockODE(x,y,e)pix = pi*x;dydx = [ y(2,:)-x/e.*y(2,:) - pi^2*cos(pix) - pix/e.*sin(pix) ];function res = shockBC(ya,yb,e)res = [ ya(1)+2yb(1)];The example passes e as an additional input argument to bvp4c.sol = bvp4c(@shockODE,@shockBC,sol,options,e);15-6915Differential Equationsbvp4c then passes this argument to the functions shockODE and shockBC whenit evaluates them. See “Additional BVP Solver Arguments” on page 15-62 formore information.2 Provide Analytical Partial Derivatives.
For this problem, the solver benefits fromusing analytical partial derivatives. The code below represents the derivativesin functions shockJac and shockBCJac.function jac = shockJac(x,y,e)jac = [ 010 -x/e ];function [dBCdya,dBCdyb] = shockBCJac(ya,yb,e)dBCdya = [ 1 00 0 ];dBCdyb = [ 0 01 0 ];shockJac and shockBCJac must accept the additional argument e, becausebvp4c passes the additional argument to all the functions the user supplies.Tell bvp4c to use these functions to evaluate the partial derivatives by settingthe options FJacobian and BCJacobian.
Also set 'Vectorized' to 'on' toindicate that the differential equation function shockODE is vectorized.options = bvpset('FJacobian',@shockJac,...'BCJacobian',@shockBCJac,...'Vectorized','on');3 Create an Initial Guess. You must provide bvp4c with a guess structure thatcontains an initial mesh and a guess for values of the solution at the meshpoints. A constant guess of y ( x ) ≡ 1 and y′ ( x ) ≡ 0 , and a mesh of five equally–2spaced points on [-1 1] suffice to solve the problem for ε = 10 . Use bvpinitto form the guess structure.sol = bvpinit([-1 -0.5 0 0.5 1],[1 0]);4 Use Continuation to Solve the Problem. To obtain the solution for the parameter–4ε = 10 , the example uses continuation by solving a sequence of problems for–2–3–4ε = 10 , 10 , 10 .
The solver bvp4c does not perform continuationautomatically, but the code's user interface has been designed to make15-70Boundary Value Problems for ODEscontinuation easy. The code uses the output sol that bvp4c produces for onevalue of e as the guess in the next iteration.e = 0.1;for i=2:4e = e/10;sol = bvp4c(@shockODE,@shockBC,sol,options,e);end5 View the Results.
Complete the example by displaying the final solutionplot(sol.x,sol.y(1,:))axis([-1 1 -2.2 2.2])title(['There is a shock at x = 0 when \epsilon = '...sprintf('%.e',e) '.'])xlabel('x')ylabel('solution y')There is a shock at x = 0 when ε =1e−004.21.51solution y0.50−0.5−1−1.5−2−1−0.8−0.6−0.4−0.20x0.20.40.60.8115-7115Differential EquationsExample: Using Continuation to Verify a Solution’s Consistent BehaviorFalkner-Skan BVPs arise from similarity solutions of viscous, incompressible,laminar flow over a flat plate. An example is2f′′′ + ff′′ + β ( 1 – ( f′ ) ) = 0for β = 0.5 on the interval [0, ∞) with boundary conditions f(0) = 0 ,f′(0) = 0 , and f′(∞) = 1 .The BVP cannot be solved on an infinite interval, and it would be impracticalto solve it for even a very large finite interval.
So, the example tries to solve asequence of problems posed on increasingly larger intervals to verify thesolution’s consistent behavior as the boundary approaches ∞ .The example imposes the infinite boundary condition at a finite point calledinfinity. The example then uses continuation in this end point to getconvergence for increasingly larger values of infinity. It uses bvpinit toextrapolate the solution sol for one value of infinity as an initial guess for thenew value of infinity. The plot of each successive solution is superimposedover those of previous solutions so they can easily be compared for consistency.Note The demo fsbvp contains the complete code for this example. The demouses subfunctions to place all required functions in a single M-file.
To run thisexample type fsbvp at the command line. See “BVP Solver Basic Syntax” onpage 15-61 and “Representing BVP Problems” on page 15-63 for moreinformation.1 Code the ODE and Boundary Condition Functions. Code the differential equationand the boundary conditions as functions that bvp4c can use.function dfdeta = fsode(eta,f)beta = 0.5;dfdeta = [ f(2)f(3)-f(1)*f(3) - beta*(1 - f(2)^2) ];function res = fsbc(f0,finf)res = [f0(1)f0(2)15-72Boundary Value Problems for ODEsfinf(2) - 1];2 Create an Initial Guess. You must provide bvp4c with a guess structure thatcontains an initial mesh and a guess for values of the solution at the meshpoints. A crude mesh of five points and a constant guess that satisfies theboundary conditions are good enough to get convergence when infinity = 3.infinity = 3;maxinfinity = 6;solinit = bvpinit(linspace(0,infinity,5),[0 0 1]);4 Solve on the Initial Interval.
The example obtains the solution for infinity = 3.It then prints the computed value of f′′(0) for comparison with the valuereported by Cebeci and Keller [2].sol = bvp4c(@fsode,@fsbc,solinit);eta = sol.x;f = sol.y;fprintf('\n');fprintf('Cebeci & Keller report that f''''(0) = 0.92768.\n')fprintf('Value computed using infinity = %g is '...'%7.5f.\n',Bnew,f(3,1))The example printsCebeci & Keller report that f''(0) = 0.92768.Value computed using infinity = 3 is 0.92915.5 Setup the Figure and Plot the Initial Solution.figureplot(eta,f(2,:),eta(end),f(2,end),'o');axis([0 maxinfinity 0 1.4]);title('Falkner-Skan equation, positive wall shear, \beta = 0.5.')xlabel('\eta')ylabel('df/d\eta')hold ondrawnowshg15-7315Differential EquationsFalkner−Skan equation, positive wall shear, β = 0.5.1.21df/dη0.80.60.40.200123η4566 Use Continuation to Solve the Problem and Plot Susbsequent Solutions.
The examplethen solves the problem for infinity = 4, 5, 6. It uses bvpinit to extrapolatethe solution sol for one value of infinity as an initial guess for the next valueof infinity. For each iteration, the example prints the computed value of f′′(0)and superimposes a plot of the solution in the existing figure.for Bnew = infinity+1:maxinfinitysolinit = bvpinit(sol,[0 Bnew]); % Extend solution to Bnew.sol = bvp4c(@fsode,@fsbc,solinit);eta = sol.x;f = sol.y;fprintf('Value computed using infinity = %g is '...'%7.5f.\n',Bnew,f(3,1))plot(eta,f(2,:),eta(end),f(2,end),'o');drawnowendhold off15-74Boundary Value Problems for ODEsThe example printsValue computed using infinity = 4 is 0.92774.Value computed using infinity = 5 is 0.92770.Value computed using infinity = 6 is 0.92770.Note that the values approach 0.92768 as reported by Cebeci and Keller.
Thesuperimposed plots confirm the consistency of the solution’s behavior.Falkner−Skan equation, positive wall shear, β = 0.5.1.21df/dη0.80.60.40.200123η456Improving BVP Solver PerformanceThe default integration properties in the BVP solver bvp4c are selected tohandle common problems. In some cases, you can improve solver performanceby changing these defaults. To do this, supply bvp4c with one or more propertyvalues in an options structure.sol = bvp4c(odefun,bcfun,solinit,options)This section:• Explains how to create, modify, and query an options structure• Describes the properties that you can use in an options structure15-7515Differential EquationsIn this and subsequent property tables, the most commonly used propertycategories are listed first, followed by more advanced categories.BVP Property CategoriesProperties CategoryProperty NamesError controlRelTol, AbsTolVectorizationVectorizedAnalytical partial derivativesFJacobian, BCJacobianMesh sizeNMaxOutput displayedStatsNote For other ways to improve solver efficiency, check “Using Continuationto Make a Good Initial Guess” on page 8-68 and the tutorial, “SolvingBoundary Value Problems for Ordinary Differential Equations in MATLABwith bvp4c,” available at ftp://ftp.mathworks.com/pub/doc/papers/bvp/.Creating and Maintaining a BVP Options StructureThe bvpset function creates an options structure that you can supply tobvp4c.
You can use bvpget to query the options structure for the value of aspecific property.Creating an Options Structure. The bvpset function accepts property name/property value pairs using the syntaxoptions = bvpset('name1',value1,'name2',value2,...)This creates a structure options in which the named properties have thespecified values. Unspecified properties retain their default values. For allproperties, it is sufficient to type only the leading characters that uniquelyidentify the property name. bvpset ignores case for property names.With no arguments, bvpset displays all property names and their possiblevalues, indicating defaults with braces {}.15-76Boundary Value Problems for ODEsModifying an Existing Options Structure. To modify an existing options argument,useoptions = bvpset(oldopts,'name1',value1,...)This overwrites any values in oldopts that are specified using name/valuepairs.
The modified structure is returned as the output argument. In the sameway, the commandoptions = bvpset(oldopts,newopts)combines the structures oldopts and newopts. In options, any values set innewopts overwrite those in oldopts.Querying an Options Structure. The bvpget function extracts a property valuefrom an options structure created with bvpset.o = bvpget(options,'name')This returns the value of the specified property, or an empty matrix [] if theproperty value is unspecified in the options structure.As with bvpset, it is sufficient to type only the leading characters that uniquelyidentify the property name; case is ignored for property names.Error Tolerance PropertiesBecause bvp4c uses a collocation formula, the numerical solution is based on amesh of points at which the collocation equations are satisfied.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
