📄 difeqdem.html
字号:
<html xmlns:mwsh="http://www.mathworks.com/namespace/mcode/v1/syntaxhighlight.dtd"> <head> <meta http-equiv="Content-Type" content="text/html; charset=utf-8"> <!--This HTML is auto-generated from an M-file.To make changes, update the M-file and republish this document. --> <title>Example: Solving a Nonlinear ODE with a Boundary Layer</title> <meta name="generator" content="MATLAB 7.1"> <meta name="date" content="2005-07-27"> <meta name="m-file" content="difeqdem"> <link rel="stylesheet" type="text/css" href="../../matlab/demos/private/style.css"> </head> <body> <div class="header"> <div class="left"><a href="matlab:edit difeqdem">Open difeqdem.m in the Editor</a></div> <div class="right"><a href="matlab:echodemo difeqdem">Run in the Command Window</a></div> </div> <div class="content"> <h1>Example: Solving a Nonlinear ODE with a Boundary Layer</h1> <introduction> <p>Illustration of toolbox use in a nontrivial problem.</p> </introduction> <h2>Contents</h2> <div> <ul> <li><a href="#1">The problem</a></li> <li><a href="#2">The approximation space</a></li> <li><a href="#3">The approximation space (continued)</a></li> <li><a href="#4">Discretization</a></li> <li><a href="#5">The numerical problem</a></li> <li><a href="#6">Linearization</a></li> <li><a href="#7">Linear system to be solved</a></li> <li><a href="#8">Linear system to be solved (continued)</a></li> <li><a href="#9">Need initial guess for y</a></li> <li><a href="#10">Iteration</a></li> <li><a href="#11">Getting ready for a smaller epsilon</a></li> <li><a href="#12">Collocation sites for new breaks</a></li> <li><a href="#13">Initial guess</a></li> <li><a href="#14">Iteration with smaller epsilon</a></li> <li><a href="#15">Very small epsilon</a></li> <li><a href="#16">Plot the breaks used for smallest epsilon</a></li> <li><a href="#17">Plot residual for smallest epsilon</a></li> </ul> </div> <h2>The problem<a name="1"></a></h2> <p>We consider the nonlinear singularly perturbed problem</p><pre> epsilon D^2g(x) + (g(x))^2 = 1 on [0..1] Dg(0) = g(1) = 0 .</pre><p>This problem is already quite difficult for epsilon = .001, so we choose a modest</p><pre class="codeinput">epsilon = .1;</pre><h2>The approximation space<a name="2"></a></h2> <p>We seek an approximate solution by collocation from C^1 piecewise cubics with a specified break sequence BREAKS, hence want the order k to be 4 and obtain the corresponding knot sequence as knots = augknt(breaks,4,2) </p><pre class="codeinput">breaks = (0:4)/4; k = 4;knots = augknt(breaks,k,2)</pre><pre class="codeoutput">knots = Columns 1 through 13 0 0 0 0 0.2500 0.2500 0.5000 0.5000 0.7500 0.7500 1.0000 1.0000 1.0000 Column 14 1.0000</pre><h2>The approximation space (continued)<a name="3"></a></h2> <p>Whatever the choice of order and knots, the corresponding spline space has dimension</p><pre class="codeinput">n = length(knots) - k</pre><pre class="codeoutput">n = 10</pre><h2>Discretization<a name="4"></a></h2> <p>The number 10 of degrees of freedom fits nicely with the fact that we expect to collocate at two sites per polynomial piece, for a total of 8 conditions, bringing us to 10 conditions altogether once we add the two side conditions. </p> <p>We choose two Gauss sites for each interval. For the `standard' interval [-1 .. 1]/2 of unit length, these are the two sites gauss = .5773502692*[-1;1]/2; From this, we obtain the whole collection of collocation sites by </p><pre class="codeinput">gauss = .5773502692*[-1;1]/2;ninterv = length(breaks)-1;temp = (breaks(2:ninterv+1)+breaks(1:ninterv))/2;temp = temp([1 1],:) + gauss*diff(breaks);colsites = temp(:).';</pre><h2>The numerical problem<a name="5"></a></h2> <p>The numerical problem we want to solve is to find a pp y of the given order and with the given knots that satisfies the (nonlinear) system </p><pre> Dy(0) = 0 (y(x))^2 + epsilon D^2y(x) = 1 for x in COLSITES y(1) = 0</pre><h2>Linearization<a name="6"></a></h2> <p>If y is our current approximation to the solution, then the linear problem for the better (?) solution z by Newton's method reads </p><pre> Dz(0) = 0 w_0(x)z(x) + epsilon D^2z(x) = b(x) for x in COLSITES z(1) = 0</pre><p>with w_0(x) := 2y(x), b(x) := (y(x))^2 + 1 .</p> <p>In fact, by choosing w_0(1) := 1, w_1(0) := 1 , and</p><pre> w_2(x) := epsilon, w_1(x) := 0 for x in COLSITES</pre><p>and choosing all other values of w_0, w_1, w_2, b not yet specified to be zero, we can give our system the uniform shape</p><pre>w_0(x)z(x) + w_1(x)Dz(x) + w_2(x)D^2z(x) = b(x) for x in SITES</pre><p>with</p><pre class="codeinput">sites = [0,colsites,1];</pre><h2>Linear system to be solved<a name="7"></a></h2> <p>This system converts into one for the B-spline coefficients of its solution z . For this, we need the zeroth, first, and second derivative at every x in SITES and for every relevant B-spline. These values are supplied by the toolbox command SPCOL. </p> <p>Here is the essential part of the online help for SPCOL:</p><pre class="codeinput"><span class="comment">%SPCOL B-spline collocation matrix.</span><span class="comment">%</span><span class="comment">% COLLOC = SPCOL(KNOTS,K,TAU) is the matrix</span><span class="comment">%</span><span class="comment">% [ D^m(i)B_j(TAU(i)) : i=1:length(TAU), j=1:length(KNOTS)-K ] ,</span><span class="comment">%</span><span class="comment">% with D^m(i)B_j the m(i)-fold derivative of B_j,</span><span class="comment">% B_j the j-th B-spline of order K for the knot sequence KNOTS,</span><span class="comment">% TAU a sequence of sites,</span><span class="comment">% both KNOTS and TAU are assumed to be nondecreasing, and</span><span class="comment">% m(i) is the integer #{ j<i : TAU(j) = TAU(i) }, i.e., the 'cumulative'</span><span class="comment">% multiplicity of TAU(i) in TAU.</span><span class="comment">%</span></pre><h2>Linear system to be solved (continued)<a name="8"></a></h2> <p>We use SPCOL to supply the matrix</p><pre> colmat = spcol(knots,k, brk2knt(sites,3) )</pre><p>with BRK2KNT used here to triple each entry of SITES, thus getting in COLMAT, for each x in SITES, value, first, and second derivative at x of all the relevant B-splines. </p> <p>From this, we get the collocation matrix by combining the row triple associated with x using the weights w_0(x), w_1(x), w_2(x) to get the row corresponding to x of the matrix of our linear system. </p><pre class="codeinput">colmat = spcol(knots,k, brk2knt(sites,3));</pre><h2>Need initial guess for y<a name="9"></a></h2> <p>We also need a current approximation y from our spline space. Initially, we get it by interpolating some reasonable initial guess from our pp space at SITES. For that guess, we use the parabola ()^2 - 1 which does satisfy the end conditions and lies in our spline space. We obtain its B-form by interpolation at SITES. We select the relevant interpolation matrix from the full matrix COLMAT. Here it is, in several (cautious) steps: </p><pre class="codeinput">intmat = colmat([2 1+(1:(n-2))*3,1+(n-1)*3],:);coefs = intmat\[0 colsites.*colsites-1 0].';y = spmak(knots,coefs.');<span class="comment">% We plot the result (it should be exactly ()^2-1 ), to be sure:</span>fnplt(y,<span class="string">'g'</span>), grid <span class="string">off</span>, axis(axis)title(<span class="string">'Initial guess (green) is ()^2-1'</span>)hold <span class="string">on</span></pre><img vspace="5" hspace="5" src="difeqdem_01.png"> <h2>Iteration<a name="10"></a></h2> <p>We can now complete the construction and solution of the linear system for the improved approximate solution z from our current guess y . In fact, with the initial guess y available, we now set up an iteration, to be terminated when the change z-y is less than a specified TOLERANCE. The max-norms of these changes are shown above. </p><pre class="codeinput">tolerance = 6.e-9;xlabel(<span class="string">'... and iterates; also the norm of the difference between iterates.'</span>)jc = -.2; hh(1) = text(.1,jc,<span class="string">'norm(z-y): '</span>);<span class="keyword">while</span> 1 vtau = fnval(y,colsites); weights=[0 1 0; [2*vtau.' zeros(n-2,1) repmat(epsilon,n-2,1)]; 1 0 0]; colloc = zeros(n,n); <span class="keyword">for</span> j=1:n colloc(j,:) = weights(j,:)*colmat(3*(j-1)+(1:3),:); <span class="keyword">end</span> coefs = colloc\[0 vtau.*vtau+1 0].'; z = spmak(knots,coefs.'); fnplt(z,<span class="string">'k'</span>) maxdif = max(max(abs(z.coefs-y.coefs))); jc = jc-.1; hh(end+1) = text(.1,jc,num2str(maxdif)); <span class="keyword">if</span> (maxdif<tolerance), <span class="keyword">break</span>, <span class="keyword">end</span> <span class="comment">% now reiterate</span> y = z;<span class="keyword">end</span></pre><img vspace="5" hspace="5" src="difeqdem_02.png"> <h2>Getting ready for a smaller epsilon<a name="11"></a></h2> <p>That looks like quadratic convergence, as expected from a Newton iteration.</p> <p>If we now decrease EPSILON, we create more of a boundary layer near the right endpoint, and this calls for a nonuniform mesh. We use NEWKNT to construct an appropriate (finer) mesh from the current approximation: </p><pre class="codeinput">knots = newknt(z, ninterv+1); breaks = knt2brk(knots);knots = augknt(breaks,4,2);n = length(knots)-k;</pre><h2>Collocation sites for new breaks<a name="12"></a></h2> <p>Next, we get the collocation sites corresponding to the new BREAKS, and then the new collocation matrix:</p><pre class="codeinput">delete(hh)ninterv = length(breaks)-1;temp = ((breaks(2:ninterv+1)+breaks(1:ninterv))/2);temp = temp([1 1], :) + gauss*diff(breaks);colsites = temp(:).';sites = [0,colsites,1];colmat = spcol(knots,k, brk2knt(sites,3));</pre><img vspace="5" hspace="5" src="difeqdem_03.png"> <h2>Initial guess<a name="13"></a></h2> <p>We obtain the initial guess y as the interpolant from the current spline space to the computed solution z . We plot the resulting interpolant (it should be close to our current solution) to be sure. </p><pre class="codeinput">intmat = colmat([2 1+(1:(n-2))*3,1+(n-1)*3],:);y = spmak(knots,[0 fnval(z,colsites) 0]/intmat.');fnplt(y,<span class="string">'g'</span>)title(<span class="string">'New initial guess (also green) for new value of epsilon'</span>), xlabel(<span class="string">''</span>)</pre><img vspace="5" hspace="5" src="difeqdem_04.png"> <h2>Iteration with smaller epsilon<a name="14"></a></h2> <p>Now we divide EPSILON by 3 and repeat the above iteration. Convergence is again quadratic.</p><pre class="codeinput">epsilon = epsilon/3; jc = -.2; hh = [];hh(1) = text(.1,jc,<span class="string">'norm(z-y): '</span>);<span class="keyword">while</span> 1 vtau = fnval(y,colsites); weights=[0 1 0; [2*vtau.' zeros(n-2,1) repmat(epsilon,n-2,1)]; 1 0 0]; colloc = zeros(n,n); <span class="keyword">for</span> j=1:n colloc(j,:) = weights(j,:)*colmat(3*(j-1)+(1:3),:); <span class="keyword">end</span> coefs = colloc\[0 vtau.*vtau+1 0].'; z = spmak(knots,coefs.'); fnplt(z,<span class="string">'b'</span>) maxdif = max(max(abs(z.coefs-y.coefs))); jc = jc-.1; hh(end+1)= text(.1,jc,num2str(maxdif)); <span class="keyword">if</span> (maxdif<tolerance), <span class="keyword">break</span>, <span class="keyword">end</span> <span class="comment">% now reiterate</span> y = z;<span class="keyword">end</span></pre><img vspace="5" hspace="5" src="difeqdem_05.png"> <h2>Very small epsilon<a name="15"></a></h2> <p>For a much smaller epsilon, we merely repeat these calculations, dividing epsilon by 3 each time.</p><pre class="codeinput">delete(hh);<span class="keyword">for</span> ee = 1:4 knots = newknt(z, ninterv+1); breaks = knt2brk(knots); knots = augknt(breaks,4,2); n = length(knots)-k; ninterv = length(breaks)-1; temp = ((breaks(2:ninterv+1)+breaks(1:ninterv))/2); temp = temp([1 1], :) + gauss*diff(breaks); colsites = temp(:).'; sites = [0,colsites,1]; colmat = spcol(knots,k, brk2knt(sites,3)); intmat = colmat([2 1+(1:(n-2))*3,1+(n-1)*3],:); y = spmak(knots,[0 fnval(z,colsites) 0]/intmat.'); fnplt(y,<span class="string">'g'</span>) epsilon = epsilon/3; <span class="keyword">while</span> 1 vtau = fnval(y,colsites); weights=[0 1 0; [2*vtau.' zeros(n-2,1) repmat(epsilon,n-2,1)]; 1 0 0]; colloc = zeros(n,n); <span class="keyword">for</span> j=1:n colloc(j,:) = weights(j,:)*colmat(3*(j-1)+(1:3),:); <span class="keyword">end</span> coefs = colloc\[0 vtau.*vtau+1 0].'; z = spmak(knots,coefs.'); fnplt(z,<span class="string">'b'</span>) maxdif = max(max(abs(z.coefs-y.coefs))); <span class="keyword">if</span> (maxdif<tolerance), <span class="keyword">break</span>, <span class="keyword">end</span> <span class="comment">% now reiterate</span> y = z; <span class="keyword">end</span><span class="keyword">end</span></pre><img vspace="5" hspace="5" src="difeqdem_06.png"> <h2>Plot the breaks used for smallest epsilon<a name="16"></a></h2> <p>Here is the final distribution of breaks, showing NEWKNT to have worked well in this case.</p><pre class="codeinput">breaks = fnbrk(fn2fm(z,<span class="string">'pp'</span>),<span class="string">'b'</span>);bb = repmat(breaks,3,1); cc = repmat([0;-1;NaN],1,length(breaks));plot(bb(:),cc(:),<span class="string">'r'</span>)title(<span class="string">'Initial guess (green) and iterates (blue) for epsilon = 1./3^j, j=2:5,'</span>)xlabel(<span class="string">'Also, the breaks used for smallest epsilon.'</span>)hold <span class="string">off</span></pre><img vspace="5" hspace="5" src="difeqdem_07.png"> <h2>Plot residual for smallest epsilon<a name="17"></a></h2> <p>Recall that we are solving the ODE</p><pre> epsilon D^2g(x) + (g(x))^2 = 1 on [0..1]</pre><p>As a check, we compute and plot the residual for the computed solution for the smallest epsilon. This, too, looks satisfactory.</p><pre class="codeinput">xx = linspace(0,1,201);plot(xx, 1 - epsilon*fnval(fnder(z,2),xx) - (fnval(z,xx)).^2)title(<span class="string">'Residual for the computed solution for smallest epsilon'</span>)</pre><img vspace="5" hspace="5" src="difeqdem_08.png"> <p class="footer">Copyright 1987-2005 C. de Boor and The MathWorks, Inc.<br> Published with MATLAB® 7.1<br></p> </div> <!--##### SOURCE BEGIN #####%% Example: Solving a Nonlinear ODE with a Boundary Layer%% Illustration of toolbox use in a nontrivial problem.% Copyright 1987-2005 C. de Boor and The MathWorks, Inc.% $Revision: 1.18.4.2 $%% The problem% We consider the nonlinear singularly perturbed problem%% epsilon D^2g(x) + (g(x))^2 = 1 on [0..1]% Dg(0) = g(1) = 0 .⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -