| blob | da811060c6d04aed29d1e029e1ab586bef3fe006 |
1 @menu
2 * Introduction to colnew::
3 * Functions and Variables for colnew::
4 * Examples for colnew::
5 * References for colnew::
6 @end menu
8 @node Introduction to colnew, Functions and Variables for colnew, Package colnew, Package colnew
9 @section Introduction to colnew
11 @var{colnew} solves mixed-order systems of boundary-value problems (BVPs)
12 in ordinary differential equations(ODEs). It is a Common Lisp translation
13 (via @var{f2cl}) of the Fortran routine COLNEW (@pxref{bader-ascher,, Bader&Ascher 1987}).
15 The method uses collocation at Gaussian points and interpolation using
16 basis functions. Approximate solutions are computed on a sequence
17 of automatically selected meshes until a user-specified set of tolerances
18 is satisfied. A damped Newton's method is used for the nonlinear iteration.
20 COLNEW has some advanced features:
22 @itemize
24 @item @b{Continuation:}
25 Most iterative numerical solvers require an initial guess. There are
26 parameter dependent problems that have an obvious solution for a special
27 value of the parameter
28 (for example a coupling constant vanishes). One can then slowly
29 vary the parameter, solving the problem numerically
30 using as a guess the previous solution, until one gets into a strong
31 coupling regime, hence the name continuation.
32 Sometimes an
33 artificial parameter is introduced to allow the use of continuation.
35 @item @b{Singularities:}
36 Another interesting feature of COLNEW is that one can manage singularities
37 in the differential equation by setting interior
38 "boundary points" on the singularities. This is because the differential
39 equation is only imposed in the interior of the intervals, but the solution
40 is checked for continuity and continuity of the derivative
41 on the boundary of the intervals.
43 @item @b{Stiff differential equations:}
44 COLNEW can be used for solving some stiff differential equations due to
45 the use of adaptive meshes.
47 @end itemize
49 The maxima interface to COLNEW exposes the full power and complexity
50 of the Fortran 77 implementation.
52 COLNEW is a modification of the package COLSYS (@pxref{ascher-1981a,,Ascher 1981a} and @ref{ascher-1981b,, Ascher 1981b}).
53 It incorporates a new basis
54 representation replacing B-splines, and improvements for
55 the linear and nonlinear algebraic equation solvers.
56 The package can be referenced as either COLNEW or COLSYS.
58 Many practical problems that are not in the standard form accepted by COLNEW
59 can be converted into this form. @xref{ascher-russell,, Asher&Russell 1981}.
61 @opencatbox{Categories:}
62 @category{Numerical methods}
63 @category{Differential equations}
64 @category{Share packages}
65 @category{Package colnew}
66 @closecatbox
68 @node Functions and Variables for colnew, Examples for colnew, Introduction to colnew, Package colnew
69 @section Functions and Variables for colnew
71 @anchor{colnew_expert}
72 @deffn {Function} colnew_expert @
73 @fname{colnew_expert} (@var{ncomp}, @var{m}, @var{aleft}, @var{aright}, @var{zeta}, @var{ipar}, @var{ltol}, @
74 @var{tol}, @var{fixpnt}, @var{ispace}, @var{fspace}, @var{iflag}, @
75 @var{fsub}, @var{dfsub}, @var{gsub}, @var{dgsub}, @var{guess})
77 @var{colnew_expert} solves mixed-order systems of boundary-value problems (BVPs) in ordinary
78 differential equations (ODEs) using a numerical collocation method.
80 @var{colnew_expert} returns the list [@var{iflag}, @var{fspace}, @var{ispace}].
81 @var{iflag} is an error flag. Lists @var{fspace} and @var{ispace} contain the
82 state of the solution
83 and can be: used by @mref{colnew_appsln} to calculate solution values
84 at arbitrary points in the solution domain; and passed back to @var{colnew_expert} to restart the solution process
85 with different arguments.
87 The function arguments are:
89 @table @code
91 @item ncomp
92 Number of differential equations (ncomp ≤ 20)
94 @item m
95 Integer list of length @var{ncomp}. m[j] is the order of the j-th
96 differential equation, with @math{1 ≤ m[j] ≤ 4}
97 and @math{mstar = sum(m[j]) ≤ 40}.
99 @item aleft
100 Left end of interval
102 @item aright
103 Right end of interval
105 @item zeta
106 Real list of length @var{mstar}. zeta[j] is the
107 j-th boundary or side condition point. The list zeta
108 must be ordered,
109 with zeta[j] ≤ zeta[j+1]. All side condition
110 points must be mesh points in all meshes used,
111 see description of ipar[11] and fixpnt below.
113 @item ipar
114 A integer list of length 11. The parameters in ipar are:
116 @itemize
117 @item @var{ipar[1]}
118 @itemize
119 @item 0, if the problem is linear
120 @item 1, if the problem is nonlinear
121 @end itemize
123 @item @var{ipar[2]} ( = k )@*
124 Number of collocation points per subinterval , where
125 @math{max(m[i]) ≤ k ≤ 7}.@*
126 If @var{ipar[2]}=0 then colnew sets k = max(max(m[i])+1, 5-max(m[i]))
128 @item @var{ipar[3]} ( = n )@*
129 Number of subintervals in the initial mesh.@*
130 If @var{ipar[3]} = 0 then colnew arbitrarily sets n = 5.
132 @item @var{ipar[4]} ( = ntol )@*
133 Number of solution and derivative tolerances.@*
134 Require 0 < @var{ntol} ≤ @var{mstar}.
136 @item @var{ipar[5]} ( = ndimf )@*
137 The length of list @var{fspace}. Its size provides a constraint on @var{nmax}.
138 Choose ipar[5] according to the formula
139 @math{ipar[5] ≥ nmax*nsizef}
140 where
141 @math{ nsizef = 4 + 3 * mstar + (5+kd) * kdm + (2*mstar-nrec) * 2*mstar}.
143 @item @var{ipar[6]} ( = ndimi )@*
144 The length of list @var{ispace}. Its size provides a constraint on @var{nmax}, the maximum
145 number of subintervals. Choose @var{ipar[6]} according to the formula
146 @math{ipar[6] ≥ nmax*nsizei}
147 where
148 @math{nsizei = 3 + kdm}
149 with
150 @math{kdm = kd + mstar}@*
151 @math{kd = k * ncomp}@*
152 @math{nrec = number of right end boundary conditions}.
154 @item @var{ipar[7]} ( = iprint )@*
155 output control
156 @itemize
157 @item -1, for full diagnostic printout
158 @item 0, for selected printout
159 @item 1, for no printout
160 @end itemize
162 @item @var{ipar[8]} ( = iread )@*
163 @itemize
164 @item 0, causes colnew to generate a uniform initial mesh.
165 @item 1,
166 if the initial mesh is provided by the user. it
167 is defined in fspace as follows: the mesh
168 aleft=x[1]<x[2]< ... <x[n]<x[n+1]=aright
169 will occupy fspace[1], ..., fspace[n+1]. the
170 user needs to supply only the interior mesh
171 points fspace[j] = x[j], j = 2, ..., n.
172 @item 2,
173 if the initial mesh is supplied by the user
174 as with @var{ipar[8]}=1, and in addition no adaptive
175 mesh selection is to be done.
176 @end itemize
178 @item @var{ipar[9]} ( = iguess )
179 @itemize
180 @item 0,
181 if no initial guess for the solution is provided
182 @item 1,
183 if an initial guess is provided by the user
184 in subroutine guess.
185 @item 2,
186 if an initial mesh and approximate solution
187 coefficients are provided by the user in @var{fspace}.
188 (the former and new mesh are the same).
189 @item 3,
190 if a former mesh and approximate solution
191 coefficients are provided by the user in @var{fspace},
192 and the new mesh is to be taken twice as coarse;
193 i.e.,every second point from the former mesh.
194 @item 4,
195 if in addition to a former initial mesh and
196 approximate solution coefficients, a new mesh
197 is provided in @var{fspace} as well.
198 (see description of output for further details
199 on iguess = 2, 3, and 4.)
200 @end itemize
202 @item @var{ipar[10]}
203 @itemize
204 @item 0, if the problem is regular
205 @item 1, if the first relax factor is =rstart, and the
206 nonlinear iteration does not rely on past covergence
207 (use for an extra sensitive nonlinear problem only).
208 @item 2, if we are to return immediately upon (a) two
209 successive nonconvergences, or (b) after obtaining
210 error estimate for the first time.
211 @end itemize
213 @item @var{ipar[11]} ( = nfxpnt , the dimension of fixpnt)@*
214 The number of fixed points in the mesh other than @var{aleft}
215 and @var{aright}.
216 The code requires that all side condition points
217 other than @var{aleft} and @var{aright} (see description of
218 @var{zeta}) must be included as fixed points in @var{fixpnt}.
219 @end itemize
222 @item ltol
223 A list of length @var{ntol=ipar[4]}. @var{ltol[j]=k} specifies
224 that the j-th tolerance in @var{tol} controls the error
225 in the k-th component of z(u).
227 The list @var{ltol} must be ordered with
228 @math{1 ≤ ltol[1] < ltol[2] < ... < ltol[ntol] ≤ mstar}.
230 @item tol
231 An list of length @var{ntol=ipar[4]}. @var{tol[j]} is the
232 error tolerance on the ltol[j]-th component
233 of z(u).
235 Thus, the code attempts to satisfy
236 for j=1,...,ntol on each subinterval
237 @math{abs(z(v)-z(u))[k] ≤ tol(j)*abs(z(u))[k]+tol(j)}
238 if v(x) is the approximate solution vector.
240 @item fixpnt
241 An list of length @var{ipar[11]}. It contains the points, other than
242 @var{aleft} and @var{aright}, which are to be included in every mesh.
243 All side condition points other than @var{aleft} and @var{aright}
244 (see @var{zeta}) be included as fixed points in @var{fixpnt}.
246 @item ispace
247 An integer work list of length @var{ipar[6]}.
249 @item fspace
250 A real work list of length @var{ipar[5]}.
252 @item fsub
253 @var{fsub} is a function f(x,z1,...,z[mstar]) which realizes
254 the system of ODEs.
256 It returns a list of @var{ncomp} values, one for each ODE.
257 Each value is the highest order
258 derivative in each ode in terms of of x,z1,...,z[mstar] .
260 @item dfsub
261 @var{dfsub} is a function df(x,z1,...,z[mstar]) for evaluating
262 the Jacobian of f.
264 @item gsub
265 Name of subroutine gsub(i,z1,z2,...,z[mstar]) for evaluating the i-th
266 component of the boundary value function g(z1,...,z[mstar]).
267 The independent variable x is not an argument of g. The value
268 @var{x=zeta[i]} must be substituted in advance.
271 @item dgsub
272 Name of subroutine dgsub(i,z1,...,z[mstar]) for evaluating the
273 i-th row of the Jacobian of g(z1,...,z[mstar]).
275 @item guess
276 Name of subroutine to evaluate the initial
277 approximation for (u(x)) and for dmval(u(x))= vector
278 of the mj-th derivatives of u(x).
279 This subroutine is needed only if using @var{ipar(9)} = 1.
281 @end table
283 The return value of @var{colnew_expert} is the list
284 @var{[iflag, fspace, ispace]}, where:
286 @table @code
288 @item iflag
289 The mode of return from @var{colnew_expert}.
290 @itemize
291 @item
292 = 1 for normal return
293 @item
294 = 0 if the collocation matrix is singular.
295 @item
296 = -1 if the expected no. of subintervals exceeds storage specifications.
297 @item
298 = -2 if the nonlinear iteration has not converged.
299 @item
300 = -3 if there is an input data error.
301 @end itemize
303 @item fspace
304 A list of floats of length @var{ipar[5]}.
306 @item ispace
307 A list of integers of length @var{ipar[6]}.
309 @end table
311 @mref{colnew_appsln} uses @var{fspace} and @var{ispace} to calculate solution values
312 at arbitrary points. The lists can also be used to restart the solution process
313 with modified meshes and parameters.
315 @opencatbox{Categories:}
316 @category{Numerical methods}
317 @category{Differential equations}
318 @category{Package colnew}
319 @closecatbox
321 @end deffn
324 @anchor{colnew_appsln}
325 @deffn {Function} colnew_appsln @
326 @fname{colnew_appsln} (@var{x}, @var{zlen}, @var{fspace}, @var{ispace})
328 Return a list of solution values from @mref{colnew_expert} results.
330 The function arguments are:
332 @table @code
334 @item x
335 List of x-coordinates to calculate solution.
336 @item zlen
337 @var{mstar}, the length of the solution list z
338 @item fspace
339 List @var{fspace} returned from @mrefdot{colnew_expert}
340 @item ispace
341 List @var{ispace} returned from @mrefdot{colnew_expert}
342 @end table
344 @opencatbox{Categories:}
345 @category{Numerical methods}
346 @category{Differential equations}
347 @category{Package colnew}
348 @closecatbox
350 @end deffn
353 @node Examples for colnew, References for colnew, Functions and Variables for colnew, Package colnew
354 @section Examples for colnew
356 COLNEW is best learned by example.
358 @subsection Example 1: A uniformly loaded beam of variable stiffness
360 The problem describes a uniformly loaded beam of variable stiffness, simply supported at both ends.
362 The problem from @ref{gawain-ball,, Gawain&Ball 1978} and is Example 1 from @ref{ascher-1981a,, Ascher 1981a}.
363 The maxima code is in file share/colnew/prob1.mac and a Fortran implementation
364 is in share/colnew/ex1.
366 @noindent The differential equation is
368 @center @math{(x^3 u'@w{}')'@w{}' = x^3 u'@w{}'@w{}'@w{}' + 6 x^2 u'@w{}'@w{}' + 6x u'@w{}' = 1} over @math{1 ≤ x ≤ 2}
370 @noindent with boundary conditions
372 @center @math{u(1) = 0, u'@w{}'(1) = 0, u(2) = 0, u'@w{}'(2) = 0}
374 @noindent The exact solution is
376 @center @math{u(x) = (1/4) (10 ln(2) - 3) (1-x) + (1/2) (1/x + (3+x) ln(x) - x)}
378 @noindent There is @var{nconc} = 1 differential equation of fourth order. The list of orders
379 @var{m} = [4] and @var{mstar} = sum(m[j]) = 4.
381 @noindent The unknown vector of length @var{mstar} is
383 @center @math{z(x) = [z_1(x),z_2(x),z_3(x),z_4(x)]}
385 @center @math{=[u(x),u'(x),u'@w{}'(x),u'@w{}'@w{}'(x)]}.
387 @noindent The differential equation is expressed as
389 @center @math{u'@w{}'@w{}'@w{}'(x) = F(x,z_1,z_2,z_3,z_4) = 1 - 6 x^2 z_3 - 6x z_2}
391 There are @var{mstar=4} boundary conditions. They are given by a
392 function @math{G(z_1,z_2,z_3,z_4)} that returns a list of length mstar.
393 The j-th boundary condition applies at @var{x = zeta[j]} and is satisfied
394 when @var{g[j] = 0}. We have
396 @c The {xxxxxx} set the column widths
397 @multitable {xxxxxxxxx} {xxxxxxxxx} {xxxxxxxxxx} {xxxxxxxxx}
398 @headitem j@ @ @ @tab zeta[j]@ @tab Condition@ @tab g[j]
399 @item 1
400 @tab 1.0
401 @tab @math{u=0}
402 @tab @math{z_1}
403 @item 2
404 @tab 1.0
405 @tab @math{u'@w{}'=0}
406 @tab @math{z_3}
407 @item 3
408 @tab 2.0
409 @tab @math{u=0}
410 @tab @math{z_1}
411 @item 4
412 @tab 2.0
413 @tab @math{u'@w{}'=0}
414 @tab @math{z_3}
415 @end multitable
417 giving @math{zeta = [1.0,1,0,2.0,2.0]}
418 and @math{G(z_1,z_2,z_3,z_4) = [z_1, z_3, z_1, z_3]}.
420 The Jacobians @var{df} and @var{dg} of @var{f} and @var{g} respectively
421 are determined symbolically.
423 The solution uses the default five collocation points per subinterval
424 and the first mesh contains only one subinterval.
425 The maximum error magnitude in the approximate solution is evaluated at 100 equidistant points
426 by function @mref{colnew_appsln} using the results returned by COLNEW and
427 compared to the estimates from the code.
429 @c ===beg===
430 @c load("colnew")$
431 @c
432 @c /* One differential equation of order 4 */
433 @c m : [4];
434 @c
435 @c /* Location of boundary conditions */
436 @c zeta : float([1,1,2,2]);
437 @c
438 @c /* Set up parameter array. Use defaults for all except as shown */
439 @c ipar : makelist(0,k,1,11);
440 @c /* initial mesh size */
441 @c ipar[3] : 1$
442 @c /* number of tolerances */
443 @c ipar[4] : 2$
444 @c /* size of real work array */
445 @c ipar[5] : 2000$
446 @c /* size of integer work array */
447 @c ipar[6] : 200$
448 @c
449 @c /* Two error tolerances (on u and its second derivative) */
450 @c ltol : [1, 3];
451 @c tol : [1d-7, 1d-7];
452 @c
453 @c /* Real work array */
454 @c fspace : makelist(0d0, k, 1, ipar[5])$
455 @c /* Integer work array */
456 @c ispace : makelist(0, k, 1, ipar[6])$
457 @c /* no internal fixed points */
458 @c fixpnt : []$
459 @c
460 @c /* Define the equations */
461 @c fsub(x, z1, z2, z3, z4) := [(1-6*x^2*z4-6*x*z3)/x^3];
462 @c df : jacobian(fsub(x,z1, z2, z3, z4),[z1,z2,z3,z4]);
463 @c define(dfsub(x, z1, z2, z3, z4), df);
464 @c
465 @c g(z1, z2, z3, z4) := [z1, z3, z1, z3];
466 @c gsub(i, z1, z2, z3, z4) :=
467 @c subst(['z1=z1,'z2=z2,'z3=z3,'z4=z4], g(z1, z2, z3, z4)[i]);
468 @c
469 @c dg:jacobian(g(z1, z2, z3, z4), [z1,z2,z3,z4]);
470 @c dgsub(i, z1, z2, z3, z4) :=
471 @c subst(['z1=z1,'z2=z2,'z3=z3,'z4=z4], row(dg, i)[1]);
472 @c
473 @c /* Exact solution */
474 @c exact(x) := [.25*(10.*log(2.)-3.)*(1.-x) + .5*(1./x+(3.+x)*log(x)-x),
475 @c -.25*(10.*log(2.)-3.) + .5*(-1./x/x+log(x)+(3.+x)/x-1.),
476 @c .5*(2./x**3+1./x-3./x/x),
477 @c .5*(-6./x**4-1./x/x+6./x**3)]$
478 @c
479 @c [iflag, fspace, ispace] :
480 @c colnew_expert(1, m, 1d0, 2d0, zeta, ipar, ltol, tol, fixpnt, ispace, fspace,
481 @c 0, fsub, dfsub, gsub, dgsub, dummy)$
482 @c
483 @c /* Calculate the error at 101 points using the known exact solution */
484 @c block([x : 1,
485 @c err : makelist(0d0, k, 1, 4),
486 @c j],
487 @c for j : 1 thru 101 do
488 @c block([],
489 @c zval : colnew_appsln([x], 4, fspace, ispace)[1],
490 @c u : float(exact(x)),
491 @c err : map(lambda([a,b], max(a,b)), err, abs(u-zval)),
492 @c x : x + 0.01),
493 @c print("The exact errors are:"),
494 @c printf(true, " ~{ ~11,4e~}~%", err));
495 @c ===end===
496 @example
497 @group
498 (%i1) load("colnew")$
500 @end group
501 @group
502 (%i2) /* One differential equation of order 4 */
503 m : [4];
505 (%o2) [4]
506 @end group
507 @group
508 (%i3) /* Location of boundary conditions */
509 zeta : float([1,1,2,2]);
511 (%o3) [1.0, 1.0, 2.0, 2.0]
512 @end group
513 @group
514 (%i4) /* Set up parameter array. Use defaults for all except as shown */
515 ipar : makelist(0,k,1,11);
516 (%o4) [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
517 @end group
518 @group
519 (%i5) /* initial mesh size */
520 ipar[3] : 1$
521 @end group
522 @group
523 (%i6) /* number of tolerances */
524 ipar[4] : 2$
525 @end group
526 @group
527 (%i7) /* size of real work array */
528 ipar[5] : 2000$
529 @end group
530 @group
531 (%i8) /* size of integer work array */
532 ipar[6] : 200$
534 @end group
535 @group
536 (%i9) /* Two error tolerances (on u and its second derivative) */
537 ltol : [1, 3];
538 (%o9) [1, 3]
539 @end group
540 @group
541 (%i10) tol : [1d-7, 1d-7];
543 (%o10) [1.0e-7, 1.0e-7]
544 @end group
545 @group
546 (%i11) /* Real work array */
547 fspace : makelist(0d0, k, 1, ipar[5])$
548 @end group
549 @group
550 (%i12) /* Integer work array */
551 ispace : makelist(0, k, 1, ipar[6])$
552 @end group
553 @group
554 (%i13) /* no internal fixed points */
555 fixpnt : []$
557 @end group
558 @group
559 (%i14) /* Define the equations */
560 fsub(x, z1, z2, z3, z4) := [(1-6*x^2*z4-6*x*z3)/x^3];
561 2
562 1 - 6 x z4 + (- 6) x z3
563 (%o14) fsub(x, z1, z2, z3, z4) := [------------------------]
564 3
565 x
566 @end group
567 @group
568 (%i15) df : jacobian(fsub(x,z1, z2, z3, z4),[z1,z2,z3,z4]);
569 [ 6 6 ]
570 (%o15) [ 0 0 - -- - - ]
571 [ 2 x ]
572 [ x ]
573 @end group
574 @group
575 (%i16) define(dfsub(x, z1, z2, z3, z4), df);
577 [ 6 6 ]
578 (%o16) dfsub(x, z1, z2, z3, z4) := [ 0 0 - -- - - ]
579 [ 2 x ]
580 [ x ]
581 @end group
582 @group
583 (%i17) g(z1, z2, z3, z4) := [z1, z3, z1, z3];
584 (%o17) g(z1, z2, z3, z4) := [z1, z3, z1, z3]
585 @end group
586 @group
587 (%i18) gsub(i, z1, z2, z3, z4) :=
588 subst(['z1=z1,'z2=z2,'z3=z3,'z4=z4], g(z1, z2, z3, z4)[i]);
590 (%o18) gsub(i, z1, z2, z3, z4) :=
591 subst(['z1 = z1, 'z2 = z2, 'z3 = z3, 'z4 = z4],
592 g(z1, z2, z3, z4) )
593 i
594 @end group
595 @group
596 (%i19) dg:jacobian(g(z1, z2, z3, z4), [z1,z2,z3,z4]);
597 [ 1 0 0 0 ]
598 [ ]
599 [ 0 0 1 0 ]
600 (%o19) [ ]
601 [ 1 0 0 0 ]
602 [ ]
603 [ 0 0 1 0 ]
604 @end group
605 @group
606 (%i20) dgsub(i, z1, z2, z3, z4) :=
607 subst(['z1=z1,'z2=z2,'z3=z3,'z4=z4], row(dg, i)[1]);
609 (%o20) dgsub(i, z1, z2, z3, z4) :=
610 subst(['z1 = z1, 'z2 = z2, 'z3 = z3, 'z4 = z4], row(dg, i) )
611 1
612 @end group
613 @group
614 (%i21) /* Exact solution */
615 exact(x) := [.25*(10.*log(2.)-3.)*(1.-x) + .5*(1./x+(3.+x)*log(x)-x),
616 -.25*(10.*log(2.)-3.) + .5*(-1./x/x+log(x)+(3.+x)/x-1.),
617 .5*(2./x**3+1./x-3./x/x),
618 .5*(-6./x**4-1./x/x+6./x**3)]$
620 @end group
621 @group
622 (%i22) [iflag, fspace, ispace] :
623 colnew_expert(1, m, 1d0, 2d0, zeta, ipar, ltol, tol, fixpnt, ispace, fspace,
624 0, fsub, dfsub, gsub, dgsub, dummy)$
627 VERSION *COLNEW* OF COLSYS .
630 THE MAXIMUM NUMBER OF SUBINTERVALS IS MIN ( 12 (ALLOWED FROM FSPACE), 16 (ALLOWED FROM ISPACE) )
632 THE NEW MESH (OF 1 SUBINTERVALS),
633 1.000000 2.000000
635 THE NEW MESH (OF 2 SUBINTERVALS),
636 1.000000 1.500000 2.000000
638 THE NEW MESH (OF 4 SUBINTERVALS),
639 1.000000 1.250000 1.500000 1.750000 2.000000
640 @end group
641 @group
642 (%i23) /* Calculate the error at 101 points using the known exact solution */
643 block([x : 1,
644 err : makelist(0d0, k, 1, 4),
645 j],
646 for j : 1 thru 101 do
647 block([],
648 zval : colnew_appsln([x], 4, fspace, ispace)[1],
649 u : float(exact(x)),
650 err : map(lambda([a,b], max(a,b)), err, abs(u-zval)),
651 x : x + 0.01),
652 print("The exact errors are:"),
653 printf(true, " ~@{ ~11,4e~@}~%", err));
654 The exact errors are:
655 1.7389E-10 6.2679E-9 2.1843E-7 9.5743E-6
656 (%o23) false
657 @end group
658 @end example
660 @subsection Example 2: Deformation of a spherical cap
662 These equations describe the small finite deformation of a thin shallow
663 spherical cap of constant thickness subject to a quadratically varying
664 axisymmetric external pressure distribution superimposed on a uniform
665 internal pressure distribution.
666 The problem is described in @ref{parker-wan,,Parker&Wan 1984} and is Example 2
667 from @ref{ascher-1981a,, Ascher 1981a}.
668 The maxima code is in file share/colnew/prob2.mac and a Fortran
669 implementation is in share/colnew/ex2.
671 There are two nonlinear differential equations
672 for @math{φ} and @math{ψ} over @math{0 < x < 1}.
674 @math{
675 (ε^4/μ)[φ'@w{}' + (1/x) φ' - (1/x^2) φ] + ψ (1-φ/x) - φ = - γ x (1-(1/2)x^2)
676 }
678 @math{ μ [ψ'@w{}' + (1/x) ψ' - (1/x^2)ψ] - φ(1-φ/(2x)) = 0 }
680 subject to boundary conditions
681 @math{φ = 0} and @math{x ψ' - 0.3 ψ + 0.7 x = 0} at x=0 and x=1.
684 For @math{ε = μ = 0.01}, two solutions exists. These are obtained by
685 starting the
686 nonlinear iteration from two different guesses to the solution:
687 initially with the default initial
688 guess; and secondly, with the initial conditions given by the
689 function @var{solutn}.
691 There are @var{nconc} = 2 differential equations of second order.
692 The list of orders @var{m} = [2,2] and
693 @var{mstar} = sum(m[i]) = 4.
695 The vector of unknowns of length @var{mstar}=4 is
696 @math{z(x) = [ φ(x), φ'(x), ψ(x), ψ'(x)]}.
698 The differential equation is expressed as
700 @math{[φ'@w{}'(x), ψ'@w{}'(x)]}
702 @math{=F(x,z_1,z_2,z_3,z_4)}
704 @math{=[z_1/x^2 - z_2/x + (z_1-z_3 (1-z_1/x) - γ x (1-x^2/2))/(ε^4/μ),
705 z_3/x^2 - z_4/x + z_1 (1-z_1/(2x))/μ]}
708 There are four boundary conditions given by list @math{zeta}
709 and function @math{G(z_1,z_2,z_3,z_4)}.
711 @multitable @columnfractions 0.1 0.1 0.4 0.4
712 @headitem j@ @ @ @tab zeta[j] @tab Condition @tab g[j]
713 @item 1
714 @tab 0.0
715 @tab @math{φ = 0}
716 @tab @math{z_1}
717 @item 2
718 @tab 0.0
719 @tab @math{x ψ' - 0.3 ψ + 0.7 x = 0}
720 @tab @math{z_3}
721 @item 3
722 @tab 1.0
723 @tab @math{φ = 0}
724 @tab @math{z_1}
725 @item 4
726 @tab 1.0
727 @tab @math{x ψ' - 0.3 ψ + 0.7 x = 0}
728 @tab @math{z_4 - 0.3@ z_3 + 0.7}
729 @end multitable
731 giving @math{zeta=[0.0,0.0,1.0,1.0]} and
732 @math{G(z_1,z_2,z_3,z_4)=[z_1, z_3, z_1, z_4-0.3*z_3+0.7]}
734 Note that @var{x} is not an argument of function @var{G}. The
735 value of @var{x=zeta[j]} must be substituted.
737 @c ===beg===
738 @c load("colnew")$
739 @c
740 @c /* Define constants */
741 @c gamma : 1.1;
742 @c eps : 0.01;
743 @c dmu : eps;
744 @c eps4mu : eps^4/dmu;
745 @c xt : sqrt(2*(gamma-1)/gamma);
746 @c /* Number of differential equations */
747 @c ncomp : 2;
748 @c /* Orders */
749 @c m : [2, 2];
750 @c
751 @c /* Interval ends */
752 @c aleft : 0.0;
753 @c aright : 1.0;
754 @c
755 @c /* Locations of side conditions */
756 @c zeta : float([0, 0, 1, 1])$
757 @c /* Set up parameter array. */
758 @c ipar : makelist(0,k,1,11);
759 @c /* Non-linear prob */
760 @c ipar[1] : 1;
761 @c /* 4 collocation points per subinterval */
762 @c ipar[2] : 4;
763 @c /* Initial uniform mesh of 10 subintervals */
764 @c ipar[3] : 10;
765 @c ipar[8] : 0;
766 @c /* Size of fspace, ispace */
767 @c ipar[5] : 40000;
768 @c ipar[6] : 2500;
769 @c /* No output */
770 @c ipar[7] : 1;
771 @c /* No initial approx is provided */
772 @c ipar[9] : 0;
773 @c /* Regular problem */
774 @c ipar[10] : 0;
775 @c /* No fixed points in mesh */
776 @c ipar[11] : 0;
777 @c /* Tolerances on all components */
778 @c ipar[4] : 4;
779 @c
780 @c /* Tolerances on all four components */
781 @c ltol : [1, 2, 3, 4];
782 @c tol : [1d-5, 1d-5, 1d-5, 1d-5];
783 @c
784 @c fspace : makelist(0d0, k, 1, ipar[5])$
785 @c ispace : makelist(0, k, 1, ipar[6])$
786 @c fixpnt : []$
787 @c
788 @c /* Define the equations */
789 @c fsub(x, z1, z2, z3, z4) :=
790 @c [z1/x/x - z2/x + (z1-z3*(1-z1/x) - gamma*x*(1-x*x/2))/eps4mu,
791 @c z3/x/x - z4/x + z1*(1-z1/2/x)/dmu];
792 @c
793 @c df : jacobian(fsub(x,z1, z2, z3, z4),[z1,z2,z3,z4]);
794 @c dfsub(x, z1, z2, z3, z4) :=
795 @c subst(['x=x,'z1=z1,'z2=z2,'z3=z3,'z4=z4], df);
796 @c
797 @c g(z1, z2, z3, z4) := [z1, z3, z1, z4 - 0.3*z3 + .7];
798 @c gsub(i, z1, z2, z3, z4) :=
799 @c subst(['z1=z1,'z2=z2,'z3=z3,'z4=z4], g(z1, z2, z3, z4)[i]);
800 @c
801 @c dg:jacobian(g(z1, z2, z3, z4), [z1,z2,z3,z4]);
802 @c dgsub(i, z1, z2, z3, z4) := subst(['z1=z1,'z2=z2,'z3=z3,'z4=z4], row(dg, i)[1]);
803 @c
804 @c /* Initial approximation function for second run */
805 @c solutn(x) :=
806 @c block([cons : gamma*x*(1-0.5*x*x),
807 @c dcons : gamma*(1-1.5*x*x),
808 @c d2cons : -3*gamma*x],
809 @c if is(x > xt) then
810 @c [[0, 0, -cons, -dcons],
811 @c [0, -d2cons]]
812 @c else
813 @c [[2*x, 2, -2*x + cons, -2 + dcons],
814 @c [0, d2cons]]);
815 @c
816 @c /* First run with default initial guess */
817 @c [iflag, fspace, ispace] :
818 @c colnew_expert(ncomp, m, aleft, aright, zeta, ipar, ltol, tol, fixpnt, ispace, fspace,
819 @c 0, fsub, dfsub, gsub, dgsub, dummy)$
820 @c
821 @c /* Check return status iflag, 1 = success */
822 @c iflag;
823 @c
824 @c /* Print values of solution at x = 0,.05,...,1 */
825 @c x : 0;
826 @c for j : 1 thru 21 do
827 @c block([],
828 @c zval : colnew_appsln([x], 4, fspace, ispace)[1],
829 @c printf(true, "~5,2f ~{~15,5e~}~%", x, zval),
830 @c x : x + 0.05);
831 @c
832 @c /* Second run with initial guess */
833 @c ipar[9] : 1;
834 @c [iflag, fspace, ispace] :
835 @c colnew_expert(ncomp, m, aleft, aright, zeta, ipar, ltol, tol, fixpnt, ispace, fspace,
836 @c 0, fsub, dfsub, gsub, dgsub, solutn)$
837 @c
838 @c /* Check return status iflag, 1 = success */
839 @c iflag;
840 @c
841 @c /* Print values of solution at x = 0,.05,...,1 */
842 @c x : 0;
843 @c for j : 1 thru 21 do
844 @c block([],
845 @c zval : colnew_appsln([x], 4, fspace, ispace)[1],
846 @c printf(true, "~5,2f ~{~15,5e~}~%", x, zval),
847 @c x : x + 0.05);
848 @c ===end===
849 @example
850 @group
851 (%i1) load("colnew")$
853 @end group
854 @group
855 (%i2) /* Define constants */
856 gamma : 1.1;
857 (%o2) 1.1
858 @end group
859 @group
860 (%i3) eps : 0.01;
861 (%o3) 0.01
862 @end group
863 @group
864 (%i4) dmu : eps;
865 (%o4) 0.01
866 @end group
867 @group
868 (%i5) eps4mu : eps^4/dmu;
869 (%o5) 1.0e-6
870 @end group
871 @group
872 (%i6) xt : sqrt(2*(gamma-1)/gamma);
873 (%o6) 0.42640143271122105
874 @end group
875 @group
876 (%i7) /* Number of differential equations */
877 ncomp : 2;
878 (%o7) 2
879 @end group
880 @group
881 (%i8) /* Orders */
882 m : [2, 2];
884 (%o8) [2, 2]
885 @end group
886 @group
887 (%i9) /* Interval ends */
888 aleft : 0.0;
889 (%o9) 0.0
890 @end group
891 @group
892 (%i10) aright : 1.0;
894 (%o10) 1.0
895 @end group
896 @group
897 (%i11) /* Locations of side conditions */
898 zeta : float([0, 0, 1, 1])$
899 @end group
900 @group
901 (%i12) /* Set up parameter array. */
902 ipar : makelist(0,k,1,11);
903 (%o12) [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
904 @end group
905 @group
906 (%i13) /* Non-linear prob */
907 ipar[1] : 1;
908 (%o13) 1
909 @end group
910 @group
911 (%i14) /* 4 collocation points per subinterval */
912 ipar[2] : 4;
913 (%o14) 4
914 @end group
915 @group
916 (%i15) /* Initial uniform mesh of 10 subintervals */
917 ipar[3] : 10;
918 (%o15) 10
919 @end group
920 @group
921 (%i16) ipar[8] : 0;
922 (%o16) 0
923 @end group
924 @group
925 (%i17) /* Size of fspace, ispace */
926 ipar[5] : 40000;
927 (%o17) 40000
928 @end group
929 @group
930 (%i18) ipar[6] : 2500;
931 (%o18) 2500
932 @end group
933 @group
934 (%i19) /* No output */
935 ipar[7] : 1;
936 (%o19) 1
937 @end group
938 @group
939 (%i20) /* No initial approx is provided */
940 ipar[9] : 0;
941 (%o20) 0
942 @end group
943 @group
944 (%i21) /* Regular problem */
945 ipar[10] : 0;
946 (%o21) 0
947 @end group
948 @group
949 (%i22) /* No fixed points in mesh */
950 ipar[11] : 0;
951 (%o22) 0
952 @end group
953 @group
954 (%i23) /* Tolerances on all components */
955 ipar[4] : 4;
957 (%o23) 4
958 @end group
959 @group
960 (%i24) /* Tolerances on all four components */
961 ltol : [1, 2, 3, 4];
962 (%o24) [1, 2, 3, 4]
963 @end group
964 @group
965 (%i25) tol : [1d-5, 1d-5, 1d-5, 1d-5];
967 (%o25) [1.0e-5, 1.0e-5, 1.0e-5, 1.0e-5]
968 @end group
969 (%i26) fspace : makelist(0d0, k, 1, ipar[5])$
970 (%i27) ispace : makelist(0, k, 1, ipar[6])$
971 @group
972 (%i28) fixpnt : []$
974 @end group
975 @group
976 (%i29) /* Define the equations */
977 fsub(x, z1, z2, z3, z4) :=
978 [z1/x/x - z2/x + (z1-z3*(1-z1/x) - gamma*x*(1-x*x/2))/eps4mu,
979 z3/x/x - z4/x + z1*(1-z1/2/x)/dmu];
981 (%o29) fsub(x, z1, z2, z3, z4) :=
982 z1 z1 x x
983 -- z1 - z3 (1 - --) + (- gamma) x (1 - ---)
984 x z2 x 2
985 [-- - -- + ----------------------------------------,
986 x x eps4mu
987 z1
988 --
989 z3 2
990 -- z1 (1 - --)
991 x z4 x
992 -- - -- + -----------]
993 x x dmu
994 @end group
995 @group
996 (%i30) df : jacobian(fsub(x,z1, z2, z3, z4),[z1,z2,z3,z4]);
997 (%o30)
998 [ z3 1 1 z1 ]
999 [ 1000000.0 (-- + 1) + -- - - 1000000.0 (-- - 1) 0 ]
1000 [ x 2 x x ]
1001 [ x ]
1002 [ ]
1003 [ z1 50.0 z1 1 1 ]
1004 [ 100.0 (1 - ---) - ------- 0 -- - - ]
1005 [ 2 x x 2 x ]
1006 [ x ]
1007 @end group
1008 @group
1009 (%i31) dfsub(x, z1, z2, z3, z4) :=
1010 subst(['x=x,'z1=z1,'z2=z2,'z3=z3,'z4=z4], df);
1012 (%o31) dfsub(x, z1, z2, z3, z4) :=
1013 subst(['x = x, 'z1 = z1, 'z2 = z2, 'z3 = z3, 'z4 = z4], df)
1014 @end group
1015 @group
1016 (%i32) g(z1, z2, z3, z4) := [z1, z3, z1, z4 - 0.3*z3 + .7];
1017 (%o32) g(z1, z2, z3, z4) := [z1, z3, z1, z4 - 0.3 z3 + 0.7]
1018 @end group
1019 @group
1020 (%i33) gsub(i, z1, z2, z3, z4) :=
1021 subst(['z1=z1,'z2=z2,'z3=z3,'z4=z4], g(z1, z2, z3, z4)[i]);
1023 (%o33) gsub(i, z1, z2, z3, z4) :=
1024 subst(['z1 = z1, 'z2 = z2, 'z3 = z3, 'z4 = z4],
1025 g(z1, z2, z3, z4) )
1026 i
1027 @end group
1028 @group
1029 (%i34) dg:jacobian(g(z1, z2, z3, z4), [z1,z2,z3,z4]);
1030 [ 1 0 0 0 ]
1031 [ ]
1032 [ 0 0 1 0 ]
1033 (%o34) [ ]
1034 [ 1 0 0 0 ]
1035 [ ]
1036 [ 0 0 - 0.3 1 ]
1037 @end group
1038 @group
1039 (%i35) dgsub(i, z1, z2, z3, z4) := subst(['z1=z1,'z2=z2,'z3=z3,'z4=z4], row(dg, i)[1]);
1041 (%o35) dgsub(i, z1, z2, z3, z4) :=
1042 subst(['z1 = z1, 'z2 = z2, 'z3 = z3, 'z4 = z4], row(dg, i) )
1043 1
1044 @end group
1045 @group
1046 (%i36) /* Initial approximation function for second run */
1047 solutn(x) :=
1048 block([cons : gamma*x*(1-0.5*x*x),
1049 dcons : gamma*(1-1.5*x*x),
1050 d2cons : -3*gamma*x],
1051 if is(x > xt) then
1052 [[0, 0, -cons, -dcons],
1053 [0, -d2cons]]
1054 else
1055 [[2*x, 2, -2*x + cons, -2 + dcons],
1056 [0, d2cons]]);
1058 (%o36) solutn(x) := block([cons : gamma x (1 - 0.5 x x),
1059 dcons : gamma (1 - 1.5 x x), d2cons : - 3 gamma x],
1060 if is(x > xt) then [[0, 0, - cons, - dcons], [0, - d2cons]]
1061 else [[2 x, 2, - 2 x + cons, - 2 + dcons], [0, d2cons]])
1062 @end group
1063 @group
1064 (%i37) /* First run with default initial guess */
1065 [iflag, fspace, ispace] :
1066 colnew_expert(ncomp, m, aleft, aright, zeta, ipar, ltol, tol, fixpnt, ispace, fspace,
1067 0, fsub, dfsub, gsub, dgsub, dummy)$
1069 @end group
1070 @group
1071 (%i38) /* Check return status iflag, 1 = success */
1072 iflag;
1074 (%o38) 1
1075 @end group
1076 @group
1077 (%i39) /* Print values of solution at x = 0,.05,...,1 */
1078 x : 0;
1079 (%o39) 0
1080 @end group
1081 @group
1082 (%i40) for j : 1 thru 21 do
1083 block([],
1084 zval : colnew_appsln([x], 4, fspace, ispace)[1],
1085 printf(true, "~5,2f ~@{~15,5e~@}~%", x, zval),
1086 x : x + 0.05);
1088 0.00 0.00000E+0 4.73042E-2 -3.39927E-32 -1.10497E+0
1089 0.05 2.36520E-3 4.73037E-2 -5.51761E-2 -1.10064E+0
1090 0.10 4.73037E-3 4.73030E-2 -1.09919E-1 -1.08765E+0
1091 0.15 7.09551E-3 4.73030E-2 -1.63796E-1 -1.06600E+0
1092 0.20 9.46069E-3 4.73039E-2 -2.16375E-1 -1.03569E+0
1093 0.25 1.18259E-2 4.73040E-2 -2.67221E-1 -9.96720E-1
1094 0.30 1.41911E-2 4.73020E-2 -3.15902E-1 -9.49092E-1
1095 0.35 1.65562E-2 4.72980E-2 -3.61986E-1 -8.92804E-1
1096 0.40 1.89215E-2 4.72993E-2 -4.05038E-1 -8.27857E-1
1097 0.45 2.12850E-2 4.72138E-2 -4.44627E-1 -7.54252E-1
1098 0.50 2.36370E-2 4.67629E-2 -4.80320E-1 -6.72014E-1
1099 0.55 2.59431E-2 4.51902E-2 -5.11686E-1 -5.81260E-1
1100 0.60 2.81093E-2 4.07535E-2 -5.38310E-1 -4.82374E-1
1101 0.65 2.99126E-2 2.98538E-2 -5.59805E-1 -3.76416E-1
1102 0.70 3.08743E-2 5.53985E-3 -5.75875E-1 -2.65952E-1
1103 0.75 3.00326E-2 -4.51680E-2 -5.86417E-1 -1.56670E-1
1104 0.80 2.55239E-2 -1.46617E-1 -5.91753E-1 -6.04539E-2
1105 0.85 1.37512E-2 -3.46952E-1 -5.93069E-1 -1.40102E-3
1106 0.90 -1.25155E-2 -7.52826E-1 -5.93303E-1 -2.86234E-2
1107 0.95 -6.94274E-2 -1.65084E+0 -5.99062E-1 -2.48115E-1
1108 1.00 2.64233E-14 1.19263E+2 -6.25420E-1 -8.87626E-1
1109 (%o40) done
1110 @end group
1111 @group
1112 (%i41) /* Second run with initial guess */
1113 ipar[9] : 1;
1114 (%o41) 1
1115 @end group
1116 @group
1117 (%i42) [iflag, fspace, ispace] :
1118 colnew_expert(ncomp, m, aleft, aright, zeta, ipar, ltol, tol, fixpnt, ispace, fspace,
1119 0, fsub, dfsub, gsub, dgsub, solutn)$
1121 @end group
1122 @group
1123 (%i43) /* Check return status iflag, 1 = success */
1124 iflag;
1126 (%o43) 1
1127 @end group
1128 @group
1129 (%i44) /* Print values of solution at x = 0,.05,...,1 */
1130 x : 0;
1131 (%o44) 0
1132 @end group
1133 @group
1134 (%i45) for j : 1 thru 21 do
1135 block([],
1136 zval : colnew_appsln([x], 4, fspace, ispace)[1],
1137 printf(true, "~5,2f ~@{~15,5e~@}~%", x, zval),
1138 x : x + 0.05);
1139 0.00 0.00000E+0 2.04139E+0 0.00000E+0 -9.03975E-1
1140 0.05 1.02070E-1 2.04139E+0 -4.52648E-2 -9.07936E-1
1141 0.10 2.04139E-1 2.04139E+0 -9.09256E-2 -9.19819E-1
1142 0.15 3.06209E-1 2.04140E+0 -1.37379E-1 -9.39624E-1
1143 0.20 4.08279E-1 2.04141E+0 -1.85020E-1 -9.67352E-1
1144 0.25 5.10351E-1 2.04152E+0 -2.34246E-1 -1.00301E+0
1145 0.30 6.12448E-1 2.04303E+0 -2.85454E-1 -1.04663E+0
1146 0.35 7.15276E-1 2.10661E+0 -3.39053E-1 -1.09916E+0
1147 0.40 8.32131E-1 -2.45181E-1 -3.96124E-1 -1.20544E+0
1148 0.45 1.77510E-2 -3.57554E+0 -4.45400E-1 -7.13543E-1
1149 0.50 2.25122E-2 1.23608E-1 -4.80360E-1 -6.70074E-1
1150 0.55 2.58693E-2 4.85257E-2 -5.11692E-1 -5.81075E-1
1151 0.60 2.80994E-2 4.11112E-2 -5.38311E-1 -4.82343E-1
1152 0.65 2.99107E-2 2.99116E-2 -5.59805E-1 -3.76409E-1
1153 0.70 3.08739E-2 5.55200E-3 -5.75875E-1 -2.65950E-1
1154 0.75 3.00325E-2 -4.51649E-2 -5.86417E-1 -1.56669E-1
1155 0.80 2.55239E-2 -1.46616E-1 -5.91753E-1 -6.04538E-2
1156 0.85 1.37512E-2 -3.46952E-1 -5.93069E-1 -1.40094E-3
1157 0.90 -1.25155E-2 -7.52826E-1 -5.93303E-1 -2.86233E-2
1158 0.95 -6.94274E-2 -1.65084E+0 -5.99062E-1 -2.48115E-1
1159 1.00 2.65413E-14 1.19263E+2 -6.25420E-1 -8.87626E-1
1160 (%o45) done
1161 @end group
1162 @end example
1164 Columns 1 (x) and 2 (@math{φ}) of the two sets of results
1165 @ifnotinfo
1166 above, and the figure below,
1167 @end ifnotinfo
1168 can be compared with Figure 1 in @ref{ascher-1981a,, Ascher 1981a}.
1170 @ifnotinfo
1171 @image{figures/colnew-ex2,8cm}
1172 @end ifnotinfo
1175 @subsection Example 3: Rotating flow of viscous incompressible fluid
1177 Example 3 from @ref{ascher-1981a,, Ascher 1981a} describes the velocities in the
1178 boundary layer produced by the rotating flow of a viscous incompressible
1179 fluid over a stationary infinite disk (@pxref{gawain-ball,,Gawain&Ball 1978}).
1181 The solution uses a number of techniques to obtain convergence.
1182 Refer to @ref{ascher-1981a,,Ascher 1981a} for details.
1184 The code is in directory share/colnew. The maxima code is in file
1185 prob3.mac. The reference Fortran implementation is in directory ex3.
1188 @subsection Example 4: Quantum Neumann equation
1190 A more sophisticated example is @ref{bellon-talon,, Bellon&Talon 2005},
1191 which deals with singularities in the
1192 solution domain, provides an initial quess to the solution
1193 and uses continuation to solve the system of non-linear
1194 differential equations.
1196 The code is in directory share/colnew. The maxima code is in file
1197 prob4.mac. The Fortran
1198 implementation is in directory ex4.
1201 @subsection Example 5: Simple example of continuation
1203 This example (@pxref{ascher-et-al,,Ascher et al@comma{} 1995@comma{} Example 9.2}) solves a numerically
1204 difficult boundary value problem using continuation.
1206 @noindent The linear differential equation is
1207 @center @math{ε u'@w{}' + x u' = -ε π^2 cos(πx) - (πx) sin(πx)}, @math{-1 < x < 1}
1209 @noindent with boundary conditions
1210 @center @math{u(-1)=-2} and @math{u(1)=0}
1212 @noindent The exact solution is
1213 @center @math{u(x) = cos(πx) + erf(x/sqrt(2ε))/erf(1/sqrt(2ε))}
1215 When @math{ε} is small the solution has a rapid transition near @math{x=0}
1216 and is difficult to solve numerically. COLNEW is able to solve the
1217 problem for directly for @math{ε=1.0e-6}, but here we will use
1218 continuation to solve it succesively for
1219 @math{ε=[1e-2,1e-3,1e-4,1e-5,1e-6]}.
1221 There is @var{nconc} = 1 differential equation of second order.
1222 The list of orders
1223 @var{m} = [2] and @var{mstar} = sum(m[j]) = 2.
1225 The unknown vector of length @var{mstar} is
1226 @math{z(x) = [z_1(x),z_2(x)] = [u(x),u'(x)]}.
1228 The differential equation is expressed as
1229 @math{[u'@w{}'(x)] = F(x,z_1,z_2) = [-(x/ε)z_2 - π^2cos(πx) - (πx/ε)sin(πx)]}
1231 There are @var{mstar=2} boundary conditions. They are given by a
1232 function @math{G(z_1,z_2)} that returns a list of length mstar.
1233 The j-th boundary condition applies at @var{x = zeta[j]} and is satisfied
1234 when @var{g[j] = 0}. We have
1236 @multitable {xxxxxxxxx} {xxxxxxxxx} {xxxxxxxxxx} {xxxxxxxxx}
1237 @headitem j@ @ @ @tab zeta[j]@ @tab Condition@ @tab g[j]
1238 @item 1
1239 @tab -1.0
1240 @tab @math{u=-2}
1241 @tab @math{z_1+2}
1242 @item 2
1243 @tab 1.0
1244 @tab @math{u=0}
1245 @tab @math{z_1}
1246 @end multitable
1248 giving @math{zeta = [-1.0,1,0]}
1249 and @math{G(z_1,z_2) = [z_1+2, z_1]}.
1251 The Jacobians @var{df} and @var{dg} of @var{f} and @var{g} respectively
1252 are determined symbolically.
1254 The ODE will be solved for multiple values of @math{ε}. The functions
1255 @var{fsub}, @var{dfsub}, @var{gsub} and @var{dgsub} are defined
1256 before @var{e} is set, so that it can be changed in the program.
1259 @c ===beg===
1260 @c load("colnew")$
1261 @c kill(e,x,z1,z2)$
1262 @c /* Exact solution */
1263 @c exact(x):=cos(%pi*x)+erf(x/sqrt(2*e))/erf(1/sqrt(2*e))$
1264 @c /* Define the equations. Do this before e is defined */
1265 @c f: [-(x/e)*z2 - %pi^2*cos(%pi*x) - (%pi*x/e)*sin(%pi*x)];
1266 @c define(fsub(x,z1,z2),f);
1267 @c df: jacobian(f,[z1,z2]);
1268 @c define(dfsub(x,z1,z2),df);
1269 @c /* Build the functions gsub and dgsub
1270 @c Use define and buildq to remove dependence on g and dg */
1271 @c g: [z1+2,z1];
1272 @c define(gsub(i,z1,z2),buildq([g],g[i]));
1273 @c dg: jacobian(g,[z1,z2]);
1274 @c define(
1275 @c dgsub(i,z1,z2),
1276 @c buildq([val:makelist(dg[i],i,1,length(dg))],block([dg:val],dg[i])));
1277 @c /* Define constant epsilon */
1278 @c e : 0.01$
1279 @c /* Number of differential equations */
1280 @c ncomp : 1$
1281 @c /* Orders */
1282 @c m : [2]$
1283 @c /* Interval ends */
1284 @c aleft:-1.0$
1285 @c aright:1.0$
1286 @c /* Locations of side conditions */
1287 @c zeta : float([-1, 1])$
1288 @c /* Set up parameter array. */
1289 @c ipar : makelist(0,k,1,11)$
1290 @c /* linear prob */
1291 @c ipar[1] : 0$
1292 @c /* 5 collocation points per subinterval */
1293 @c ipar[2] : 5$
1294 @c /* Initial uniform mesh of 1 subintervals */
1295 @c ipar[3] : 1$
1296 @c ipar[8] : 0$
1297 @c /* Size of fspace, ispace */
1298 @c ipar[5] : 10000$
1299 @c ipar[6] : 1000$
1300 @c /* No output. Don't do this for development. */
1301 @c ipar[7]:1$
1302 @c /* No initial guess is provided */
1303 @c ipar[9] : 0$
1304 @c /* Regular problem */
1305 @c ipar[10] : 0$
1306 @c /* No fixed points in mesh */
1307 @c ipar[11] : 0$
1308 @c /* Tolerances on two components */
1309 @c ipar[4] : 2$
1310 @c /* Two error tolerances (on u and its derivative)
1311 @c Relatively large tolerances to keep the example small */
1312 @c ltol : [1, 2]$
1313 @c tol : [1e-4, 1e-4]$
1314 @c fspace : makelist(0e0, k, 1, ipar[5])$
1315 @c ispace : makelist(0, k, 1, ipar[6])$
1316 @c fixpnt : []$
1317 @c /* First run with default initial guess.
1318 @c Returns iflag. 1 = success */
1319 @c ([iflag, fspace, ispace] :
1320 @c colnew_expert(ncomp, m, aleft, aright, zeta, ipar, ltol, tol,
1321 @c fixpnt, ispace, fspace,
1322 @c 0, fsub, dfsub, gsub, dgsub, dummy),
1323 @c if (iflag#1) then error("On return from colnew_expert: iflag = ",iflag),
1324 @c iflag);
1325 @c /* Function to generate equally spaced list of values */
1326 @c xlist(xmin,xmax,n):=block([dx:(xmax-xmin)/n],makelist(i,i,0,n)*dx+xmin)$
1327 @c /* x values for solution. Cluster around x=0 */
1328 @c X: xlist(aleft,aright,500)$
1329 @c /* Generate solution values for z1=u(x) */
1330 @c ans:colnew_appsln(X,2,fspace,ispace)$
1331 @c z:maplist(first,ans)$
1332 @c Z:[z]$
1333 @c /* Compare with exact solution and report */
1334 @c y:float(map(exact,X))$
1335 @c maxerror:apply(max,abs(y-z));
1336 @c printf(true," e: ~8,3e iflag ~3d Mesh size ~3d max error ~8,3e~%",
1337 @c e,iflag,ispace[1],maxerror);
1338 @c /* Now use continuation to solve for progressively smaller e
1339 @c Use previous solution as initial guess and every second point
1340 @c from previous mesh as initial mesh */
1341 @c ipar[9] : 3$
1342 @c /* Run COLNEW using continuation for new value of e
1343 @c Set new mesh size ipar[3] from previous size ispace[1]
1344 @c Push list of values of z1=u(x) on to list Z */
1345 @c run_it(e_):=block(
1346 @c e:e_,
1347 @c ipar[3]:ispace[1],
1348 @c [iflag, fspace, ispace]:
1349 @c colnew_expert(ncomp,m,aleft,aright,zeta,ipar,ltol,tol,fixpnt,
1350 @c ispace,fspace,0,fsub,dfsub,gsub,dgsub,dummy),
1351 @c if (iflag#1) then error("On return from colnew_expert: iflag =",iflag),
1352 @c ans:colnew_appsln(X,2,fspace,ispace),
1353 @c z:maplist(first,ans),
1354 @c push(z,Z),
1355 @c y:float(map(exact,X)),
1356 @c maxerror:apply(max,abs(y-z)),
1357 @c printf(true," e: ~8,3e iflag ~3d Mesh size ~3d max error ~8,3e~%",
1358 @c e,iflag,ispace[1],maxerror),
1359 @c iflag
1360 @c )$
1361 @c for e_ in [1e-3,1e-4,1e-5,1e-6] do run_it(e_)$
1362 @c /* Z is list of solutions z1 = u(x). Restore order. */
1363 @c Z:reverse(Z)$
1364 @c /* Plot z1=u(x) for each value of e
1365 @c plot2d([
1366 @c [discrete,X,Z[1]], [discrete,X,Z[2]], [discrete,X,Z[3]],
1367 @c [discrete,X,Z[4]], [discrete,X,Z[5]]],
1368 @c [legend,"e=1e-2","e=1e-3","e=1e-4","e=1e-5","e=1e-6"],
1369 @c [xlabel,"x"],[ylabel,"u(x)"],
1370 @c [png_file,"./colnew-ex5.png"]); */
1371 @c done$
1372 @c ===end===
1373 @example
1374 (%i1) load("colnew")$
1375 (%i2) kill(e,x,z1,z2)$
1376 @group
1377 (%i3) /* Exact solution */
1378 exact(x):=cos(%pi*x)+erf(x/sqrt(2*e))/erf(1/sqrt(2*e))$
1379 @end group
1380 @group
1381 (%i4) /* Define the equations. Do this before e is defined */
1382 f: [-(x/e)*z2 - %pi^2*cos(%pi*x) - (%pi*x/e)*sin(%pi*x)];
1383 x z2 %pi x sin(%pi x) 2
1384 (%o4) [- ---- - ---------------- - %pi cos(%pi x)]
1385 e e
1386 @end group
1387 @group
1388 (%i5) define(fsub(x,z1,z2),f);
1389 x z2 %pi x sin(%pi x)
1390 (%o5) fsub(x, z1, z2) := [- ---- - ----------------
1391 e e
1392 2
1393 - %pi cos(%pi x)]
1394 @end group
1395 @group
1396 (%i6) df: jacobian(f,[z1,z2]);
1397 [ x ]
1398 (%o6) [ 0 - - ]
1399 [ e ]
1400 @end group
1401 @group
1402 (%i7) define(dfsub(x,z1,z2),df);
1403 [ x ]
1404 (%o7) dfsub(x, z1, z2) := [ 0 - - ]
1405 [ e ]
1406 @end group
1407 @group
1408 (%i8) /* Build the functions gsub and dgsub
1409 Use define and buildq to remove dependence on g and dg */
1410 g: [z1+2,z1];
1411 (%o8) [z1 + 2, z1]
1412 @end group
1413 @group
1414 (%i9) define(gsub(i,z1,z2),buildq([g],g[i]));
1415 (%o9) gsub(i, z1, z2) := [z1 + 2, z1]
1416 i
1417 @end group
1418 @group
1419 (%i10) dg: jacobian(g,[z1,z2]);
1420 [ 1 0 ]
1421 (%o10) [ ]
1422 [ 1 0 ]
1423 @end group
1424 @group
1425 (%i11) define(
1426 dgsub(i,z1,z2),
1427 buildq([val:makelist(dg[i],i,1,length(dg))],block([dg:val],dg[i])));
1428 (%o11) dgsub(i, z1, z2) := block([dg : [[1, 0], [1, 0]]], dg )
1429 i
1430 @end group
1431 @group
1432 (%i12) /* Define constant epsilon */
1433 e : 0.01$
1434 @end group
1435 @group
1436 (%i13) /* Number of differential equations */
1437 ncomp : 1$
1438 @end group
1439 @group
1440 (%i14) /* Orders */
1441 m : [2]$
1442 @end group
1443 @group
1444 (%i15) /* Interval ends */
1445 aleft:-1.0$
1446 @end group
1447 (%i16) aright:1.0$
1448 @group
1449 (%i17) /* Locations of side conditions */
1450 zeta : float([-1, 1])$
1451 @end group
1452 @group
1453 (%i18) /* Set up parameter array. */
1454 ipar : makelist(0,k,1,11)$
1455 @end group
1456 @group
1457 (%i19) /* linear prob */
1458 ipar[1] : 0$
1459 @end group
1460 @group
1461 (%i20) /* 5 collocation points per subinterval */
1462 ipar[2] : 5$
1463 @end group
1464 @group
1465 (%i21) /* Initial uniform mesh of 1 subintervals */
1466 ipar[3] : 1$
1467 @end group
1468 (%i22) ipar[8] : 0$
1469 @group
1470 (%i23) /* Size of fspace, ispace */
1471 ipar[5] : 10000$
1472 @end group
1473 (%i24) ipar[6] : 1000$
1474 @group
1475 (%i25) /* No output. Don't do this for development. */
1476 ipar[7]:1$
1477 @end group
1478 @group
1479 (%i26) /* No initial guess is provided */
1480 ipar[9] : 0$
1481 @end group
1482 @group
1483 (%i27) /* Regular problem */
1484 ipar[10] : 0$
1485 @end group
1486 @group
1487 (%i28) /* No fixed points in mesh */
1488 ipar[11] : 0$
1489 @end group
1490 @group
1491 (%i29) /* Tolerances on two components */
1492 ipar[4] : 2$
1493 @end group
1494 @group
1495 (%i30) /* Two error tolerances (on u and its derivative)
1496 Relatively large tolerances to keep the example small */
1497 ltol : [1, 2]$
1498 @end group
1499 (%i31) tol : [1e-4, 1e-4]$
1500 (%i32) fspace : makelist(0e0, k, 1, ipar[5])$
1501 (%i33) ispace : makelist(0, k, 1, ipar[6])$
1502 (%i34) fixpnt : []$
1503 @group
1504 (%i35) /* First run with default initial guess.
1505 Returns iflag. 1 = success */
1506 ([iflag, fspace, ispace] :
1507 colnew_expert(ncomp, m, aleft, aright, zeta, ipar, ltol, tol,
1508 fixpnt, ispace, fspace,
1509 0, fsub, dfsub, gsub, dgsub, dummy),
1510 if (iflag#1) then error("On return from colnew_expert: iflag = ",iflag),
1511 iflag);
1512 (%o35) 1
1513 @end group
1514 @group
1515 (%i36) /* Function to generate equally spaced list of values */
1516 xlist(xmin,xmax,n):=block([dx:(xmax-xmin)/n],makelist(i,i,0,n)*dx+xmin)$
1517 @end group
1518 @group
1519 (%i37) /* x values for solution. Cluster around x=0 */
1520 X: xlist(aleft,aright,500)$
1521 @end group
1522 @group
1523 (%i38) /* Generate solution values for z1=u(x) */
1524 ans:colnew_appsln(X,2,fspace,ispace)$
1525 @end group
1526 (%i39) z:maplist(first,ans)$
1527 (%i40) Z:[z]$
1528 @group
1529 (%i41) /* Compare with exact solution and report */
1530 y:float(map(exact,X))$
1531 @end group
1532 @group
1533 (%i42) maxerror:apply(max,abs(y-z));
1534 (%o42) 6.881499912125832e-7
1535 @end group
1536 @group
1537 (%i43) printf(true," e: ~8,3e iflag ~3d Mesh size ~3d max error ~8,3e~%",
1538 e,iflag,ispace[1],maxerror);
1539 e: 1.000E-2 iflag 1 Mesh size 16 max error 6.881E-7
1540 (%o43) false
1541 @end group
1542 @group
1543 (%i44) /* Now use continuation to solve for progressively smaller e
1544 Use previous solution as initial guess and every second point
1545 from previous mesh as initial mesh */
1546 ipar[9] : 3$
1547 @end group
1548 @group
1549 (%i45) /* Run COLNEW using continuation for new value of e
1550 Set new mesh size ipar[3] from previous size ispace[1]
1551 Push list of values of z1=u(x) on to list Z */
1552 run_it(e_):=block(
1553 e:e_,
1554 ipar[3]:ispace[1],
1555 [iflag, fspace, ispace]:
1556 colnew_expert(ncomp,m,aleft,aright,zeta,ipar,ltol,tol,fixpnt,
1557 ispace,fspace,0,fsub,dfsub,gsub,dgsub,dummy),
1558 if (iflag#1) then error("On return from colnew_expert: iflag =",iflag),
1559 ans:colnew_appsln(X,2,fspace,ispace),
1560 z:maplist(first,ans),
1561 push(z,Z),
1562 y:float(map(exact,X)),
1563 maxerror:apply(max,abs(y-z)),
1564 printf(true," e: ~8,3e iflag ~3d Mesh size ~3d max error ~8,3e~%",
1565 e,iflag,ispace[1],maxerror),
1566 iflag
1567 )$
1568 @end group
1569 @group
1570 (%i46) for e_ in [1e-3,1e-4,1e-5,1e-6] do run_it(e_)$
1571 e: 1.000E-3 iflag 1 Mesh size 20 max error 3.217E-7
1572 e: 1.000E-4 iflag 1 Mesh size 40 max error 3.835E-7
1573 e: 1.000E-5 iflag 1 Mesh size 38 max error 8.690E-9
1574 e: 1.000E-6 iflag 1 Mesh size 60 max error 6.313E-7
1575 @end group
1576 @group
1577 (%i47) /* Z is list of solutions z1 = u(x). Restore order. */
1578 Z:reverse(Z)$
1579 @end group
1580 @group
1581 (%i48) /* Plot z1=u(x) for each value of e
1582 plot2d([
1583 [discrete,X,Z[1]], [discrete,X,Z[2]], [discrete,X,Z[3]],
1584 [discrete,X,Z[4]], [discrete,X,Z[5]]],
1585 [legend,"e=1e-2","e=1e-3","e=1e-4","e=1e-5","e=1e-6"],
1586 [xlabel,"x"],[ylabel,"u(x)"],
1587 [png_file,"./colnew-ex5.png"]); */
1588 done$
1589 @end group
1590 @end example
1592 @ifnotinfo
1593 The figure below shows the solution for
1594 @math{ε=[10^{-2},10^{-3},10^{-4},10^{-5},10^{-6}]}.
1596 @image{figures/colnew-ex5,8cm}
1597 @end ifnotinfo
1602 @node References for colnew, , Examples for colnew, Package colnew
1603 @section References for colnew
1605 @itemize
1607 @item @anchor{gawain-ball}
1608 (Gawain&Ball 1978) T. H. Gawain and R. E. Ball,
1609 Improved Finite Difference Formulas for Boundary Value Problems,
1610 International Journal for Numerical Methods in Engineering 12, no. 7 (1978)
1611 1151–60.
1612 @url{https://doi.org/10.1002/nme.1620120706, doi:10.1002/nme.1620120706}
1615 @item @anchor{ascher-1979a}
1616 (Ascher 1979a) U. Ascher, J. Christiansen and R. D. Russell,
1617 A collocation solver for mixed order systems of boundary value problems,
1618 Math. Comp. 33 (1979), 659-679,
1619 @url{https:/doi.org/10.1090/S0025-5718-1979-0521281-7,
1620 doi:10.1090/S0025-5718-1979-0521281-7}
1622 @item @anchor{ascher-1979b}
1623 (Ascher 1979b) U. Ascher, J. Christiansen and R. D. Russell,
1624 COLSYS - a collocation code for boundary value problems,
1625 in @i{Codes for boundary-value problems in ordinary differential equations},
1626 Lecture Notes in Computer Science 76, Springer Verlag,
1627 B. Childs et. al. (eds.) (1979), 164-185,
1628 ISBN 978-3-540-09554-5
1630 @item @anchor{ascher-1981a}
1631 (Ascher 1981a) U. Ascher, J. Christiansen and R. D. Russell,
1632 Collocation software for boundary-value odes,
1633 ACM Trans. Math Software 7 (1981), 209-222.
1634 @url{https:/doi.org/10.1145/355945.355950,
1635 doi:10.1145/355945.355950}
1637 @item @anchor{ascher-1981b}
1638 (Ascher 1981b) U. Ascher, U., J. Christiansen, and R. D. Russell.
1639 ‘Algorithm 569: COLSYS: Collocation Software for Boundary-Value ODEs [D2]’.
1640 ACM Transactions on Mathematical Software 7, no. 2 (June 1981): 223–29.
1641 @url{https:/doi.org/10.1145/355945.355951,
1642 doi:10.1145/355945.355951}
1644 @item @anchor{ascher-russell}
1645 (Ascher&Russell 1981) U. Ascher and R. D. Russell.
1646 ‘Reformulation of Boundary Value Problems into “Standard” Form’.
1647 SIAM Review 23, no. 2 (April 1981), 238–54,
1648 @url{https:/doi.org/10.1137/1023039, doi:10.1137/1023039}
1650 @item @anchor{parker-wan}
1651 (Parker&Wan 1984) David F. Parker and Frederic Y. M. Wan,
1652 ‘Finite Polar Dimpling of Shallow Caps Under Sub-Buckling Axisymmetric
1653 Pressure Distributions’.
1654 SIAM Journal on Applied Mathematics 44, no. 2 (April 1984): 301–26,
1655 @url{https://doi.org/10.1137/0144022, doi:10.1137/0144022}
1657 @item @anchor{bader-ascher}
1658 (Bader&Ascher 1987) G. Bader and U. Ascher,
1659 A new basis implementation for a mixed order boundary value ode solver,
1660 SIAM J. Scient. Stat. Comput. (1987), 483-500,
1661 @url{https://doi.org/10.1137/0908047, doi:10.1137/0908047}
1663 @item @anchor{ascher-et-al}
1664 (Ascher et al, 1995) Uri M. Ascher, Robert M. M. Mattheij,
1665 and Robert D. Russell.
1666 Numerical Solution of Boundary Value Problems for Ordinary Differential
1667 Equations.
1668 Classics in Applied Mathematics 13, SIAM, (1995),
1669 ISBN 978-0-89871-354-1
1671 @item @anchor{bellon-talon}
1672 (Bellon&Talon 2005) M. Bellon, M. Talon,
1673 Spectrum of the quantum Neumann model,
1674 Physics Letters A, Volume 337, Issues 4–6, pp 360-368,
1675 @url{https://doi.org/10.1016/j.physleta.2005.02.002,
1676 doi:10.1016/j.physleta.2005.02.002}
1677 @url{https://arxiv.org/abs/hep-th/0407005,arXiv:hep-th/0407005}
1678 @end itemize