8 years, 40 days

## What is the true equation your time and...

Maple

This post is closely related to this one Centered Divided Difference approximations whose purpose is to build Finite-Difference (FD) approxmation schemes.
In this latter post I also talked, without explicitely naming it, about Truncation Error, see for instance https://en.wikiversity.org/wiki/Numerical_Analysis/Truncation_Errors.

I am foccusing here on a less known concept named "Equivalent Equation" (sometimes named "Modified Equation").
The seminal paper (no free acccess, which is surprising since it was first published 50 years ago) is this one by Warming and Hyett https://www.sciencedirect.com/science/article/pii/0021999174900114.
For a scholar example you can see here https://guillod.org/teaching/m2-b004/TD2-solution.pdf.
An alternative method to that of Warming and Hyett to derive the Equivalent Equation is given here (in French)
http://www.numdam.org/item/M2AN_1997__31_4_459_0.pdf (Carpentier et al.)

I never heard of the concept of Modified Equation applied to elliptic PDEs ; it's main domain of application is advection PDEs (parabolic PDEs in simpler cases).

Basically, any numerical scheme for solving an ODE or a PDE has a Truncation Error: this means there is some discrepancy between the exact solution of this PDE and the solution the scheme provides.
This discrepancy depends on the truncation error, in space alone for ODEs or elliptic PDEs, or in space and time for hyperbolic or advection PDEs.

One main problem with the Truncation Error is that it doesn't enable to understand the impact it will have on the solution it returns. For instance, will this sheme introduce diffusion, antidiffusion, scattering, ...
The aim of the Modified Equation is to answer these questions by constructing the continuous equation a given numerical scheme solves with a zero truncation error.
This is the original point of view Warming and Hyett developped in their 1974 paper.

This is a subject I work on 30 years ago (Maple V), and later in 2010.
It is very easy with Maple to do the painstaking development that Warming and Hyett did by hand half a century ago. And it is even so easy that the trick used by Carpentier et al. to make the calculations easier has lost some of its interest.

Two examples are given in the attched file
EquivalentEquation.mw

If a moderator thinks that this post should be merged with Centered Divided Difference approximations, he can do it freely, I won't be offended.

## Centered Divided Difference approximatio...

Maple

This code enables building Centered Divided Difference (CDD) approximations of derivatives of a univariate function.
Depending on the stencil we choose we can get arbitrary high order approximations.

The extension to bivariate functions is based upon what is often named tensorization in numerical analysis: for instance diff(f(x, y), [x, y] is obtained this way (the description here is purely notional)

1. Let CDD_x the CDD approximation of diff( f(x), x) ) .
CDD_x is a linear combination of shifted replicates of f(x)
2. Let s one of this shifted replicates
Let CDD_y(s) the CDD approximation of diff( s(y), y) ) .
3. Replace in CDD_x all shifted replicates by their corresponding expression CDD_y(s)

REMARKS:

• When I write for instance "approximation of diff(f(x), x)", this must be intended as a short for "approximation of diff(f(x), x) at point x=a"
• I agree that a notation like, for instance, diff(f(a), a) is not rigourous and that something like a Liebnitz notation would be better. Unfortunately I don't know how to get it when I use mtaylor.

 > restart:

CDDF stands for Cendered Divided Difference Formula

 > CDDF := proc(f, A, H, n, stencil)   description "f = target function,\nA = point where the derivatives are approximated,\nH = step,\nn = order of the derivative,\nstencil = list of points for the divided differenceCDDF\n";   local tay, p, T, sol, unknown, Unknown, expr:   tay := (s, m) -> convert(                      eval(                        convert(                          taylor(op(0, f)(op(1, f)), op(1, f)=A, m),                          Diff                        ),                        op(1, f)=A+s*H),                      polynom                    ):   p   := numelems(stencil):   T   := add(alpha[i]*tay(i, p+1), i in stencil):   T   := convert(%, diff):   if p > n+1 then     sol := solve([seq(coeff(T, h, i)=0, i in subsop(n+1=NULL, [$0..p]))], [seq(alpha[i], i in stencil)])[]; else sol := solve([seq(coeff(T, H, i)=0, i in subsop(n+1=NULL, [$0..n]))], [seq(alpha[i], i in stencil)])[];   end if:   if and(is~(rhs~(sol)=~0)[]) then     WARNING("no solution found"):     return   else     unknown := union(indets~(rhs~(sol))[])[];     Unknown := simplify(solve(eval(T, sol) = Diff(op(0, f)(A), A$n), unknown)); sol := lhs~(sol) =~ eval(rhs~(sol), unknown=Unknown); expr := normal(eval(add(alpha[i]*op(0, f)(A+i*H), i in stencil), sol)); end if: return expr end proc:  > Describe(CDDF)  # f = target function, # A = point where the derivatives are approximated, # H = # step, # n = order of the derivative, # stencil = list of points for the divided # differenceCDDF # CDDF( f, A, H, n, stencil )  > # 2-point approximation of diff(f(x), x) | x=a CDDF(f(x), a, h, 1, [-1, 1]); convert(simplify(mtaylor(%, h=0, 4)), Diff);  (1)  > # 3-point approximation of diff(f(x), x$2) | x=a CDDF(f(x), a, h, 2, [-1, 0, 1]); convert(simplify(mtaylor(%, h=0)), Diff);
 (2)
 > # 5-point pproximation of diff(f(x), x$2) | x=a CDDF(f(x), a, h, 2, [$-2..2]); convert(simplify(mtaylor(%, h=0, 8)), Diff);
 (3)
 > # 7-point approximation of diff(f(x), x$2) | x=a CDDF(f(x), a, h, 2, [$-3..3]); # simplify(taylor(%, h=0, 10)); convert(simplify(mtaylor(%, h=0, 10)), Diff);
 (4)
 > # 4-point staggered approximation of diff(f(x), x$3) | x=a CDDF(f(x), a, h, 3, [seq(-3/2..3/2, 1)]); convert(simplify(mtaylor(%, h=0, 6)), Diff);  (5)  > # 6-point staggered approximation of diff(f(x), x$3) | x=a CDDF(f(x), a, h, 3, [seq(-5/2..5/2, 1)]); # simplify(taylor(%, h=0, 8)); convert(simplify(mtaylor(%, h=0, 8)), Diff);
 (6)
 > # 5-point approximation of diff(f(x), x$4) | x=a CDDF(f(x), a, h, 4, [$-2..2]); convert(simplify(mtaylor(%, h=0, 8)), Diff);
 (7)
 > # 7-point approximation of diff(f(x), x$4) | x=a CDDF(f(x), a, h, 4, [$-3..3]); convert(simplify(mtaylor(%, h=0, 10)), Diff);
 (8)

A FEW 2D EXTENSIONS

 > # diff(f(x, y), [x, y]) approximation over a (2 by 2)-point stencil stencil := [-1, 1]: # step 1: approximate diff(f(x, y), x) over stencil < stencil > fx  := CDDF(f(x), a, h, 1, stencil): fx  := eval(% , f=(u -> f[u](y))): ix  := [indets(fx, function)[]]: # step 2: approximate diff(g(y), y) over stencil < stencil > where #         g represents any function in fx. fxy := add(map(u -> CDDF(u, b, k, 1, stencil)*coeff(fx, u), ix)): # step 3: rewrite fxy in a more convenient form [seq(u=op([0, 0], u)(op([0, 1], u), op(1, u)), u in indets(fxy, function))]: fxy := simplify( eval(fxy, %) ); convert(mtaylor(fxy, [h=0, k=0]), Diff)
 (9)
 > # Approximation of diff(f(x, y), [x, x, y, y] a (3 by 3)-point stencil stencil := [-1, 0, 1]: # step 1: approximate diff(f(x, y), x) over stencil < stencil > fx  := CDDF(f(x), a, h, 2, stencil): fx  := eval(% , f=(u -> f[u](y))): ix  := [indets(fx, function)[]]: # step 2: approximate diff(g(y), y) over stencil < stencil > where #         g represents any function in fx. fxy := add(map(u -> CDDF(u, b, k, 2, stencil)*coeff(fx, u), ix)): # step 3: rewrite fxy in a more convenient form [seq(u=op([0, 0], u)(op([0, 1], u), op(1, u)), u in indets(fxy, function))]: fxy := simplify( eval(fxy, %) ); convert(mtaylor(fxy, [h=0, k=0], 8), Diff)
 (10)
 > # Approximation of diff(f(x, y), [x, x, y] a (3 by 2)-point stencil stencil_x := [-1, 0, 1]: stencil_y := [-1, 1]: # step 1: approximate diff(f(x, y), x) over stencil < stencil > fx  := CDDF(f(x), a, h, 2, stencil_x): fx  := eval(% , f=(u -> f[u](y))): ix  := [indets(fx, function)[]]: # step 2: approximate diff(g(y), y) over stencil < stencil > where #         g represents any function in fx. fxy := add(map(u -> CDDF(u, b, k, 1, stencil_y)*coeff(fx, u), ix)): # step 3: rewrite fxy in a more convenient form [seq(u=op([0, 0], u)(op([0, 1], u), op(1, u)), u in indets(fxy, function))]: fxy := simplify( eval(fxy, %) ); convert(mtaylor(fxy, [h=0, k=0], 6), Diff)
 (11)
 > # Approximation of the laplacian of f(x, y) stencil := [-1, 0, 1]: # step 1: approximate diff(f(x, y), x) over stencil < stencil > fx  := CDDF(f(x), a, h, 2, stencil): fy  := CDDF(f(y), b, k, 2, stencil): fxy := simplify( eval(fx, f=(u -> f(u, b))) + eval(fy, f=(u -> f(a, u))) ); convert(mtaylor(fxy, [h=0, k=0], 6), Diff)
 (12)
 >

## Maple proof of Poncelet's theorem for th...

Maple

This post is inspired by a serie of questions from @JAMET.
I wondered if it was possible to prove plane geometry theorems with the geometry package.

Here is an illustration for the Poncelet's theorem for the triangle (French designation), see for instance
https://en.wikipedia.org/wiki/Poncelet%27s_closure_theorem

Are any of you interested in challenging the geometry package to proof other plane geometry theorems?

 > restart:

Poncelet's theorem for the triangle

Let ABC a triangle, c its incircle (center omega, radius r) and C its circumcircle (center Omega, radius R).
Let D the distance between omega and Omega.

then R^2 - D^2 - 2*r*R = 0

Proof

Without loss of generality one sets :

x(A) = y(A) = 0
and  y(B) = 0

ABC is a non degerated triangle provided x(B) <> 0 and y(C) <> 0

 > with(geometry):
 > kernelopts(version);
 (1)
 > assume(x__B <> 0): assume(y__C <> 0):
 > point(A, 0, 0); point(B, x__B, 0); point(C, x__C, y__C);
 (2)
 > triangle(T, [A, B, C])
 (3)
 > bisector(bA, A, T); bisector(bB, B, T); eA := isolate(Equation(bA, [x, y]), y): eB := isolate(Equation(bB, [x, y]), y): xc := solve(rhs(eA)=rhs(eB), x): yc := eval(rhs(eA), x=xc): point(omega, xc, yc): r := distance(omega, line(lAB, [A, B]))
 (4)
 > circumcircle(C, T, 'centername' = Omega); R := radius(C);
 (5)
 > Oo := distance(Omega, omega)
 (6)
 > S := simplify(R^2 - Oo^2 - 2*r*R)  assuming x__B::real, x__C::real, y__C::real
 (7)
 > simplify(S) assuming x__B > 0, y__C > 0; simplify(S) assuming x__B > 0, y__C < 0; simplify(S) assuming x__B < 0, y__C > 0; simplify(S) assuming x__B < 0, y__C < 0;
 (8)
 >

Improvements of the geometry package

It already appears that (some) assumptions are not (always) correctly taken into account. This is a weak point which requires, as in the attached mw, to use an indirect approache to construct the incircle.

As a matter of fact, the procedure incircle, whose first lines are

showstat(incircle)

geometry:-incircle := proc(inci, T)
local cname, d, A, B, l1, l2, dis, x, y, tmp, msg;
1   if nargs < 2 or 3 < nargs then
2     error "wrong number of arguments"
end if;
3   if geometry:-form(T) <> ('triangle2d') then
4     error "wrong type of arguments"
end if;
5   if nargs = 3 and op(1,args[3]) = ('centername') and type(op(2,args[3]),'name') then
6     cname := op(2,args[3])
else
7     cname := cat('center_',inci)
end if;
8   if geometry:-method(T) = (':-points') then
9     d := geometry:-DefinedAs(T);
10     A := op(1,d);
11     B := op(2,d);
12     msg := sprintf("find the bisector of %a at vertex %a",T,A);
13     userinfo(2,geometry,msg);
14     geometry:-bisector('l1',A,T);
15     msg := sprintf("find the bisector of %a at vertex %a",T,B);
16     userinfo(2,geometry,msg);
17     geometry:-bisector('l2',B,T);
18     msg := sprintf("find the intersection of the two bisectors");
19     userinfo(2,geometry,msg);
20     geometry:-intersection(cname,l1,l2);


requires that the two bissectors are not parallel(call to geometry:-intersection).

Since the non-parallelism of bisectors is trivial for all non-degenerate triangles, why doesn't incircle inherit this property rather than not being able to decide if the bisectors are parallel or not?)

Here is a detail of what happens and the endless loop in which incircle seems to be caught

kernelopts(version);
Maple 2015.2, APPLE UNIVERSAL OSX, Dec 20 2015, Build ID 1097895

AreParallel(bA, bB, 'condition'):
AreParallel: hint: cannot determine if -y__C*(x__B^2)^(1/2)*(x__C*(x__B^2)^(1/2)-x__B*(x__B^2)^(1/2)-x__B*((x__B-x__C)^2+y__C^2)^(1/2))+y__C*(x__B^2)^(1/2)*(x__B*(x__C^2+y__C^2)^(1/2)+x__C*(x__B^2)^(1/2)) is zero

assume(lhs(condition) <> 0);
AreParallel(bA, bB);
false
intersection(J, bA, bB);
J

infolevel[geometry] := 4:
incircle(inc, T);
incircle: find the bisector of T at vertex A
incircle: find the bisector of T at vertex B
incircle: find the intersection of the two bisectors
intersection: find the intersection between two lines l1 and l2
intersection: one point of intersection
incircle: find the radius of the incircle
line: define the line from two points
circle: define the circle from its center and radius
circle: hint: abs(x__B^2*y__C/(csgn(x__B)*x__B+(x__B^2-2*x__B*x__C+x__C^2+y__C^2)^(1/2)+(x__C^2+y__C^2)^(1/2)))/(x__B^2)^(1/2) > 0
Error, (in geometry:-circle) not enough information: the radius might not be positive
assume(abs(x__B^2*y__C/(csgn(x__B)*x__B+(x__B^2-2*x__B*x__C+x__C^2+y__C^2)^(1/2)+(x__C^2+y__C^2)^(1/2)))/(x__B^2)^(1/2) > 0):

assume(abs(x__B^2*y__C/(csgn(x__B)*x__B+(x__B^2-2*x__B*x__C+x__C^2+y__C^2)^(1/2)+(x__C^2+y__C^2)^(1/2)))/(x__B^2)^(1/2) > 0):
incircle(inc, T);
incircle: find the bisector of T at vertex A
incircle: find the bisector of T at vertex B
incircle: find the intersection of the two bisectors
intersection: find the intersection between two lines l1 and l2
AreParallel: hint: cannot determine if -y__C*(x__B^2)^(1/2)*(x__C*(x__B^2)^(1/2)-x__B*(x__B^2)^(1/2)-x__B*((x__B-x__C)^2+y__C^2)^(1/2))+y__C*(x__B^2)^(1/2)*(x__B*(x__C^2+y__C^2)^(1/2)+x__C*(x__B^2)^(1/2)) is zero
intersection: two given lines intersect each other if -y__C*(x__B^2)^(1/2)*(x__C*(x__B^2)^(1/2)-x__B*(x__B^2)^(1/2)-x__B*((x__B-x__C)^2+y__C^2)^(1/2))+y__C*(x__B^2)^(1/2)*(x__B*(x__C^2+y__C^2)^(1/2)+x__C*(x__B^2)^(1/2)) <> 0
Error, (in geometry:-intersection) not enough information

assume(-y__C*(x__B^2)^(1/2)*(x__C*(x__B^2)^(1/2)-x__B*(x__B^2)^(1/2)-x__B*((x__B-x__C)^2+y__C^2)^(1/2))+y__C*(x__B^2)^(1/2)*(x__B*(x__C^2+y__C^2)^(1/2)+x__C*(x__B^2)^(1/2)) <> 0):
incircle(inc, T);
incircle: find the bisector of T at vertex A
incircle: find the bisector of T at vertex B
incircle: find the intersection of the two bisectors
intersection: find the intersection between two lines l1 and l2
intersection: one point of intersection
incircle: find the radius of the incircle
line: define the line from two points
circle: define the circle from its center and radius
circle: hint: abs(x__B^2*y__C/(csgn(x__B)*x__B+(x__C^2+y__C^2)^(1/2)+(x__B^2-2*x__B*x__C+x__C^2+y__C^2)^(1/2)))/(x__B^2)^(1/2) > 0
Error, (in geometry:-circle) not enough information: the radius might not be positive



## Suggestion: improvement of help pages in...

Maple

A simple suggestion...

I would appreciate being able to open multiple help pages simultaneously instead of just one.
This seems to me particularly interesting when you have to browse back and forth between several related items.

## A Maple generator of LatTeX \tabular str...

Maple

In two recent questions

dnaviaux  raised concerns about saving informations he could include in a typeset report.

In the first question these informations were mainly printf terminal outputs, while the second question was oriented to informations of type Table (in the sense of DocumentTools:-Layout:-Table).

Acer has given a direct answer to this last question.
On my side I proposed a less direct one based on a LaTeXapproach.

Since several years I use to generate a LaTeX source code for typeset reports directly from Maple (this is not very difficult once you know LaTeX).
These reports often contain data tables Mbut latex(M)is not the good way to code these data tables into LaTeX.
That is why I turned towards the LaTeX/tabularstructure

\begin{tabular}{...}
...
\end{tabular]

Data tables are very simple structures and converting them programmatically from Maple Matrix into  LaTeX tabular is rather simple: the most important thing you must remember is that some characters have a special role in LaTeX  (for instance "\" , "_" and "$") while having a special meaning in Maple too: so beware of conflicts. The procedure Tabular in the attached file takes as argument a Matrix and returns the LaTeX code of its tabular counterpart. Note that it could be easily extendend to accept a DataFrameas input. This is a very simple simple version and fonts (\rm, \bf, \it... or others), styles (\scriptsize, \large, ...) and much more customization features could be accounted for without any difficulty. For this particular example inspired by dnaviaux first question, it would be probably better to replace the LaTeX code by one of thes two: \begin{tabular}{|c @{$\pm$} c|} .... Asdfsdg & +5.190e+01 & 2.5950e+00 \\  or : \begin{tabular}{|c|c|} .... Asdfsdg & \multicolumn{2}{c|}{$+8.680e-01 \pm 0.9190} \\  dnaviaux second question concerned a DocumentTools:-Layout:-Table object. For what I understood this table has a simple structure (no rowspan neither columnspan used). But, knowing that tabular can manage merged columns (\multicolumn) and merged rows (\multirow, provided the ad hoc pasckage is used), I wondered if it would be possible to generate the LaTeX/tabular code corresponding to a Table using rowspan and columnspan? This is done by the procedure Tabular_Table in the attached file. Here again this is a very simple procedure which doesn't exploit all the informations a Table structure contains (backgroundstylefillcoloralignalignmentseparator [Tabular_Table separates all columns and rows] ...). The "raw" rendering is not bad but could be improved: • adjust the way \multirowcenters the content of a cell (Xis badly placed in example 3) • adjust the spaces between line/clineand the text • by using the package cellspace • or my modufying \arraystretch with \renewcommand. • one could also add a legend with \caption • parameterize Tabular_Table to accept other column separators • manage the colors (for instance a blue for a math expression and a black for a text) • ... Finally Tabularand Tabular_Tablecould have an optional argument file: {file::{symbol, string}:= terminal}  set by default to terminalor wich could be the name of a .texfile This work is still under developpement and I would be happy to exchange with you on this topic. Happy New Year to all of you Tabular_from_Maple.mw Here is the content of the .tex file (I used Welcome-to-CoCalc.texto create/compile it) \documentclass{article} % set font encoding for PDFLaTeX, XeLaTeX, or LuaTeX \usepackage{ifxetex,ifluatex} \if\ifxetex T\else\ifluatex T\else F\fi\fi T% \usepackage{fontspec} \else \usepackage[T1]{fontenc} \usepackage[utf8]{inputenc} \usepackage{lmodern} \fi \usepackage{hyperref} % DO NOT FORGET THIS PACKAGE !!! \usepackage{multirow} \title{Tabulars from Maple} \author{mmcdara} \begin{document} \maketitle \begin{itemize} \item[\bf{Example 1}] \begin{tabular}{|c|c|c|} \hline \hline Quantity & Nominal value & Uncertainty \\ \hline Asdfsdg & +9.060e+01 & 4.5298e+00 \\ \hline Bdfg & +1.437e+01 & 7.1870e-01 \\ \hline C123 & +8.025e+01 & 4.0125e+00 \\ \hline Ddf sdfg dsfg & +9.614e+00 & 4.8072e-01 \\ \hline \hline \end{tabular} \item[]\vspace*{10mm} \item[\bf{Example 3}] \begin{tabular}{| c | c | c | c | c |} \hline \multicolumn{1}{|c|}{\multirow{4}{*}{X}} & \multicolumn{2}{c|}{\multirow{1}{*}{A1}} & \multicolumn{1}{c|}{\multirow{2}{*}{A2}} & \multicolumn{1}{c|}{\multirow{1}{*}{A3}} \\ \cline{2-2} \cline{3-3} \cline{5-5} & \multicolumn{1}{c|}{\multirow{2}{*}{B1}} &\multicolumn{1}{c|}{\multirow{1}{*}{B2}} & & \multicolumn{1}{c|}{\multirow{3}{*}{B3}} \\ \cline{3-3} \cline{4-4} & & \multicolumn{2}{c|}{{\cos \left( \omega\,t+\phi \right) }^{\mathstrut}_{\atop{}{\mathstrut}}$} & \\ \cline{2-2} \cline{3-3} \cline{4-4} & \multicolumn{3}{c|}{${{\frac {1}{\Gamma  \left( x \right) }\sqrt {{{\rm e}^{-{\frac {{t}^{2}}{\pi }}}}}}}^{\mathstrut}_{\atop{}{\mathstrut}}\$} &  \\
\hline
\end{tabular}

\end{itemize}

\end{document}

Image of the resulting PDF:

For comparison here is the Tableof Example 3 that Mapledisplays in the worksheet

 1 2 3 4 5 6 Page 2 of 6
﻿