Skip to content
Permalink

Comparing changes

Choose two branches to see what’s changed or to start a new pull request. If you need to, you can also or learn more about diff comparisons.

Open a pull request

Create a new pull request by comparing changes across two branches. If you need to, you can also . Learn more about diff comparisons here.
base repository: mockingbirdnest/Principia
Failed to load repositories. Confirm that selected base ref is valid, then try again.
Loading
base: fbf6ff8344f5
Choose a base ref
...
head repository: mockingbirdnest/Principia
Failed to load repositories. Confirm that selected head ref is valid, then try again.
Loading
compare: 069007944353
Choose a head ref
  • 3 commits
  • 5 files changed
  • 1 contributor

Commits on Sep 16, 2020

  1. Copy the full SHA
    f1640f3 View commit details
  2. merge

    eggrobin committed Sep 16, 2020
    Copy the full SHA
    3149fa2 View commit details

Commits on Sep 17, 2020

  1. Merge pull request #2720 from eggrobin/gauss-legendre

    Tables for Gauss-Legendre quadrature
    eggrobin authored Sep 17, 2020

    Verified

    This commit was created on GitHub.com and signed with GitHub’s verified signature. The key has expired.
    Copy the full SHA
    0690079 View commit details
Showing with 2,694 additions and 0 deletions.
  1. +107 −0 mathematica/gauss_legendre.wl
  2. +1,290 −0 numerics/gauss_legendre_weights.mathematica.h
  3. +1,290 −0 numerics/legendre_roots.mathematica.h
  4. +2 −0 numerics/numerics.vcxproj
  5. +5 −0 numerics/numerics.vcxproj.filters
107 changes: 107 additions & 0 deletions mathematica/gauss_legendre.wl
Original file line number Diff line number Diff line change
@@ -0,0 +1,107 @@
(* ::Package:: *)

(* ::Text:: *)
(*See https://dlmf.nist.gov/3.5#v.*)


ClearAll[LegendreNodes];
LegendreNodes[n_]:=LegendreNodes[n]=Sort[N[Last/@Flatten[Solve[LegendreP[n,x]==0,x]],300]]


MonicLegendreP[n_][x_]:=LegendreP[n,x]/Coefficient[LegendreP[n,x],x^n]


ClearAll[LegendreWeights];
LegendreWeights[n_]:=LegendreWeights[n]=Table[NIntegrate[MonicLegendreP[n][x]/((x-LegendreNodes[n][[k]])MonicLegendreP[n]'[LegendreNodes[n][[k]]]),{x,-1,1},AccuracyGoal->\[Infinity],PrecisionGoal->100,WorkingPrecision->200,Method->{"GlobalAdaptive",Method->"GaussKronrodRule"}],{k,1,n}]


nodes=Table[
Table[
{n,k,Re[N[LegendreNodes[n][[k+1]],46]]},
{k,0,n-1}],
{n,0,50}];


nodes[[;;,;;,-1]]//TableForm


And@@Flatten[Table[nodes[[n+1,k+1,3]]==-nodes[[n+1,-(k+1),3]],{n,0,50},{k,0,n-1}]]


weights=Table[
Table[
{n,k,Re[N[LegendreWeights[n][[k+1]],46]]},
{k,0,n-1}],
{n,0,50}];


weights[[;;,;;,-1]]//TableForm


And@@Flatten[Table[weights[[n+1,k+1,3]]==weights[[n+1,-(k+1),3]],{n,0,50},{k,0,n-1}]]


decimalFloatLiteral[x_Real|x:0,exponentWidth_Integer,signed_]:=
With[
{m=MantissaExponent[x][[1]]*10,
e=If[x==0,0,MantissaExponent[x][[2]]-1]},
StringJoin[
StringRiffle[#,"'"]&/@
{{If[signed&&m>=0,"+",Nothing]},{#[[1]]},
If[Length[#]>1,{"."},Nothing],
If[Length[#]>1,StringPartition[#[[2]],UpTo[5]],Nothing],
If[e!=0,"e"<>If[e>0,"+","-"]<>IntegerString[e,10,exponentWidth],Nothing]}&[
StringSplit[ToString[m],"."]]]]


SetDirectory[NotebookDirectory[]]


Export[
"..\\numerics\\legendre_roots.mathematica.h",
"
#pragma once
#include \"numerics/fixed_arrays.hpp\"
namespace principia {
namespace numerics {
// Roots of the Legendre polynomials of degree up to 50.
constexpr FixedStrictlyLowerTriangularMatrix<double, "<>ToString[51]<>">
LegendreRoots{{{
"<>Flatten@Map[
With[
{n=#[[1]],k=#[[2]],z=#[[3]]},
" /*"<>If[k==0,"n="<>StringPadLeft[ToString[n],2]<>", "," "]<>"k="<>StringPadLeft[ToString[k],2]<>"*/"<>decimalFloatLiteral[z,1,True]<>",\n"]&,
nodes,{2}]<>"}}};
} // namespace numerics
} // namespace principia
",
"text"]


Export[
"..\\numerics\\gauss_legendre_weights.mathematica.h",
"
#pragma once
#include \"numerics/fixed_arrays.hpp\"
namespace principia {
namespace numerics {
// Weights for Gauss-Legendre quadrature with up to 50 points.
constexpr FixedStrictlyLowerTriangularMatrix<double, "<>ToString[51]<>">
GaussLegendreWeights{{{
"<>Flatten@Map[
With[
{n=#[[1]],k=#[[2]],z=#[[3]]},
" /*"<>If[k==0,"n="<>StringPadLeft[ToString[n],2]<>", "," "]<>"k="<>StringPadLeft[ToString[k],2]<>"*/"<>decimalFloatLiteral[z,1,False]<>",\n"]&,
weights,{2}]<>"}}};
} // namespace numerics
} // namespace principia
",
"text"]
Loading