Age

120 Reputation

6 Badges

10 years, 38 days

MaplePrimes Activity


These are replies submitted by Age

@acer  Thanks!

This works quite well. I'll try and implement it into the rest of the code. 

 

Thanks again!

@acer thanks for the reply! 

I had trialed a few integration methods, however the function wasn't an ArrayInterpolation in these trials, it was a function given out by dsolve. Anyway, in those trials, I found the 2D integration CubaCuhre to be the best. Anyway, I have created a simplified version of what I am trying to accomplish. 

Thanks again. 

 


``

``

``

Sample Integration Sheet

 

 

 

Currently Using this to convert grid to spline

 

Instead of x and y coords I'm usng t and r. So TCut is max range in the t-direction, similarly for R. nt is the grid size in the t direction and nr is the size in the r direction.

This function just builds spline fit and follows this: https://www.maplesoft.com/support/help/maple/view.aspx?path=examples%2FInterpolation_and_Smoothing

 

 

 

GridToSpline := proc (initGrid, TCut, RCut, nt, nr) local it, ir, VectT, VectR, Lt, Lr, currentGrid, SplineGrid; VectT := Vector(nt); VectR := Vector(nr); Lt := 1/(nt-1); Lr := 1/(nr-1); for it to nt do VectT[it] := TCut*((2*it-2)*Lt-1) end do; for ir to nr do VectR[ir] := RCut*(ir-1)*Lr end do; currentGrid := Array(1 .. nt, 1 .. nr, proc (i, j) options operator, arrow; evalf(initGrid[i][j]) end proc, datatype = float[8]); SplineGrid := proc (a, b) options operator, arrow; CurveFitting:-ArrayInterpolation([VectT, VectR], Array(currentGrid), Array(1 .. 1, 1 .. 1, 1 .. 2, [[[a, b]]]), 'method' = 'spline')[1, 1] end proc; return SplineGrid end proc:
NULL

``

NULL

 

Here is an example of a grid I'm trying to solve.

 

 

TCut := .3;

TCut := .3

RCut := .3
nt := 29
nr := 29

 

 

 

SampleGrid := Matrix(29, 29, {(1, 1) = -2.69866813533281, (1, 2) = -2.69866820985674, (1, 3) = -2.69866827617016, (1, 4) = -2.69866837749853, (1, 5) = -2.69866850212862, (1, 6) = -2.69866863715926, (1, 7) = -2.69866877073077, (1, 8) = -2.69866889363839, (1, 9) = -2.69866900008058, (1, 10) = -2.69866908759731, (1, 11) = -2.69866915646738, (1, 12) = -2.69866920891275, (1, 13) = -2.69866924843909, (1, 14) = -2.69866927960774, (1, 15) = -2.69866930859281, (1, 16) = -2.69866934517168, (1, 17) = -2.69866940758526, (1, 18) = -2.69866953354284, (1, 19) = -2.69866980479120, (1, 20) = -2.69867040196861, (1, 21) = -2.69867172735934, (1, 22) = -2.69867468019801, (1, 23) = -2.69868127520480, (1, 24) = -2.69869603438140, (1, 25) = -2.69872912299730, (1, 26) = -2.69880342434324, (1, 27) = -2.69897051873145, (1, 28) = -2.69934681269525, (1, 29) = -2.70019531250000, (2, 1) = -2.69866603047084, (2, 2) = -2.69866625169965, (2, 3) = -2.69866644691331, (2, 4) = -2.69866674348742, (2, 5) = -2.69866710514694, (2, 6) = -2.69866749259564, (2, 7) = -2.69866787060424, (2, 8) = -2.69866821286631, (2, 9) = -2.69866850389706, (2, 10) = -2.69866873830742, (2, 11) = -2.69866891844754, (2, 12) = -2.69866905155222, (2, 13) = -2.69866914728661, (2, 14) = -2.69866921625670, (2, 15) = -2.69866926986058, (2, 16) = -2.69866932201265, (2, 17) = -2.69866939402006, (2, 18) = -2.69866952574703, (2, 19) = -2.69866980038930, (2, 20) = -2.69867039952334, (2, 21) = -2.69867172602142, (2, 22) = -2.69867467947623, (2, 23) = -2.69868127482055, (2, 24) = -2.69869603417946, (2, 25) = -2.69872912289265, (2, 26) = -2.69880342429014, (2, 27) = -2.69897051870583, (2, 28) = -2.69934681268504, (2, 29) = -2.70019531250000, (3, 1) = -2.69865211033863, (3, 2) = -2.69865336241077, (3, 3) = -2.69865445573850, (3, 4) = -2.69865610484201, (3, 5) = -2.69865809438942, (3, 6) = -2.69866019585749, (3, 7) = -2.69866221080965, (3, 8) = -2.69866399852749, (3, 9) = -2.69866548417810, (3, 10) = -2.69866665085402, (3, 11) = -2.69866752273294, (3, 12) = -2.69866814665308, (3, 13) = -2.69866857702353, (3, 14) = -2.69866886614376, (3, 15) = -2.69866906000736, (3, 16) = -2.69866919897384, (3, 17) = -2.69866932332963, (3, 18) = -2.69866948588351, (3, 19) = -2.69866977829260, (3, 20) = -2.69867038746753, (3, 21) = -2.69867171953946, (3, 22) = -2.69867467603810, (3, 23) = -2.69868127301993, (3, 24) = -2.69869603324795, (3, 25) = -2.69872912241713, (3, 26) = -2.69880342405215, (3, 27) = -2.69897051859236, (3, 28) = -2.69934681264019, (3, 29) = -2.70019531250000, (4, 1) = -2.69857559192165, (4, 2) = -2.69858297079499, (4, 3) = -2.69858934181473, (4, 4) = -2.69859887732956, (4, 5) = -2.69861024891079, (4, 6) = -2.69862207755826, (4, 7) = -2.69863320751553, (4, 8) = -2.69864286743489, (4, 9) = -2.69865069942731, (4, 10) = -2.69865668685300, (4, 11) = -2.69866103505642, (4, 12) = -2.69866405363369, (4, 13) = -2.69866606808988, (4, 14) = -2.69866736808164, (4, 15) = -2.69866818667951, (4, 16) = -2.69866870083484, (4, 17) = -2.69866904479461, (4, 18) = -2.69866933294639, (4, 19) = -2.69866969570361, (4, 20) = -2.69867034354260, (4, 21) = -2.69867169650278, (4, 22) = -2.69867466411125, (4, 23) = -2.69868126691862, (4, 24) = -2.69869603016244, (4, 25) = -2.69872912087587, (4, 26) = -2.69880342329623, (4, 27) = -2.69897051823817, (4, 28) = -2.69934681250191, (4, 29) = -2.70019531250000, (5, 1) = -2.69815364103587, (5, 2) = -2.69819775736826, (5, 3) = -2.69823537129639, (5, 4) = -2.69829118183365, (5, 5) = -2.69835687584068, (5, 6) = -2.69842403321791, (5, 7) = -2.69848588072185, (5, 8) = -2.69853822709065, (5, 9) = -2.69857949015361, (5, 10) = -2.69861008904377, (5, 11) = -2.69863160905534, (5, 12) = -2.69864606076830, (5, 13) = -2.69865538171212, (5, 14) = -2.69866118567008, (5, 15) = -2.69866469370480, (5, 16) = -2.69866676914122, (5, 17) = -2.69866799700861, (5, 18) = -2.69866877449387, (5, 19) = -2.69866940276336, (5, 20) = -2.69867019209140, (5, 21) = -2.69867161923136, (5, 22) = -2.69867462516171, (5, 23) = -2.69868124750408, (5, 24) = -2.69869602058751, (5, 25) = -2.69872911620663, (5, 26) = -2.69880342105693, (5, 27) = -2.69897051720914, (5, 28) = -2.69934681210573, (5, 29) = -2.70019531250000, (6, 1) = -2.69580300166118, (6, 2) = -2.69607068008258, (6, 3) = -2.69629574928807, (6, 4) = -2.69662647000063, (6, 5) = -2.69701006252744, (6, 6) = -2.69739450225433, (6, 7) = -2.69773989523952, (6, 8) = -2.69802385072690, (6, 9) = -2.69824051704782, (6, 10) = -2.69839565518402, (6, 11) = -2.69850083935603, (6, 12) = -2.69856887877845, (6, 13) = -2.69861113639229, (6, 14) = -2.69863647184190, (6, 15) = -2.69865120559918, (6, 16) = -2.69865955858620, (6, 17) = -2.69866421319792, (6, 18) = -2.69866682171521, (6, 19) = -2.69866840997105, (6, 20) = -2.69866969415176, (6, 21) = -2.69867137253520, (6, 22) = -2.69867450429699, (6, 23) = -2.69868118889116, (6, 24) = -2.69869599243597, (6, 25) = -2.69872910282171, (6, 26) = -2.69880341478773, (6, 27) = -2.69897051438702, (6, 28) = -2.69934681103519, (6, 29) = -2.70019531250000, (7, 1) = -2.68255270905730, (7, 2) = -2.68420269183551, (7, 3) = -2.68556927039621, (7, 4) = -2.68755569167303, (7, 5) = -2.68982174085213, (7, 6) = -2.69204188519929, (7, 7) = -2.69398007273551, (7, 8) = -2.69552020727773, (7, 9) = -2.69665146213758, (7, 10) = -2.69742910118661, (7, 11) = -2.69793453459984, (7, 12) = -2.69824780099437, (7, 13) = -2.69843426151571, (7, 14) = -2.69854146697746, (7, 15) = -2.69860129800890, (7, 16) = -2.69863385101707, (7, 17) = -2.69865119943290, (7, 18) = -2.69866033490219, (7, 19) = -2.69866522085018, (7, 20) = -2.69866814551416, (7, 21) = -2.69867062880846, (7, 22) = -2.69867415067918, (7, 23) = -2.69868102227875, (7, 24) = -2.69869591459725, (7, 25) = -2.69872906677682, (7, 26) = -2.69880339831608, (7, 27) = -2.69897050713020, (7, 28) = -2.69934680832502, (7, 29) = -2.70019531250000, (8, 1) = -2.60686981106159, (8, 2) = -2.61720825621594, (8, 3) = -2.62563727540733, (8, 4) = -2.63774342800024, (8, 5) = -2.65129960389903, (8, 6) = -2.66424133250227, (8, 7) = -2.67516577737618, (8, 8) = -2.68350330362916, (8, 9) = -2.68935662535048, (8, 10) = -2.69319142854833, (8, 11) = -2.69556414442440, (8, 12) = -2.69696419297891, (8, 13) = -2.69775825071771, (8, 14) = -2.69819388604856, (8, 15) = -2.69842625642532, (8, 16) = -2.69854727740723, (8, 17) = -2.69860905128715, (8, 18) = -2.69864009741297, (8, 19) = -2.69865562174525, (8, 20) = -2.69866364141685, (8, 21) = -2.69866853562029, (8, 22) = -2.69867318623206, (8, 23) = -2.69868058132789, (8, 24) = -2.69869571442844, (8, 25) = -2.69872897658248, (8, 26) = -2.69880335813610, (8, 27) = -2.69897048981909, (8, 28) = -2.69934680196416, (8, 29) = -2.70019531250000, (9, 1) = -2.16846353937698, (9, 2) = -2.23425128880023, (9, 3) = -2.28707572809903, (9, 4) = -2.36196107905338, (9, 5) = -2.44412091008287, (9, 6) = -2.52028432208616, (9, 7) = -2.58206991833073, (9, 8) = -2.62698442443030, (9, 9) = -2.65683978853650, (9, 10) = -2.67530598933600, (9, 11) = -2.68608729120433, (9, 12) = -2.69209644788596, (9, 13) = -2.69532209803656, (9, 14) = -2.69700104276086, (9, 15) = -2.69785291472636, (9, 16) = -2.69827602403542, (9, 17) = -2.69848245044229, (9, 18) = -2.69858169930828, (9, 19) = -2.69862895750863, (9, 20) = -2.69865157471866, (9, 21) = -2.69866311747724, (9, 22) = -2.69867077012410, (9, 23) = -2.69867951050463, (9, 24) = -2.69869524249446, (9, 25) = -2.69872876980332, (9, 26) = -2.69880326838836, (9, 27) = -2.69897045202630, (9, 28) = -2.69934678830801, (9, 29) = -2.70019531250000, (10, 1) = .399782631544292, (10, 2) = -0.223230825286887e-1, (10, 3) = -.357156479975148, (10, 4) = -.825507840664823, (10, 5) = -1.32860914943849, (10, 6) = -1.78026685157640, (10, 7) = -2.12983346183128, (10, 8) = -2.36923077886601, (10, 9) = -2.51800244358535, (10, 10) = -2.60380139320853, (10, 11) = -2.65054483715776, (10, 12) = -2.67492220656520, (10, 13) = -2.68721033766193, (10, 14) = -2.69324013010359, (10, 15) = -2.69613566822986, (10, 16) = -2.69750183210721, (10, 17) = -2.69813712499983, (10, 18) = -2.69842905264433, (10, 19) = -2.69856199874022, (10, 20) = -2.69862239583270, (10, 21) = -2.69865047391492, (10, 22) = -2.69866531830100, (10, 23) = -2.69867716974714, (10, 24) = -2.69869424135636, (10, 25) = -2.69872834337068, (10, 26) = -2.69880308809385, (10, 27) = -2.69897037783243, (10, 28) = -2.69934676194892, (10, 29) = -2.70019531250000, (11, 1) = 15.2509102692799, (11, 2) = 12.6264291046438, (11, 3) = 10.5462806670587, (11, 4) = 7.66265417268524, (11, 5) = 4.61175757088330, (11, 6) = 1.94893912064126, (11, 7) = -0.655913599774656e-2, (11, 8) = -1.25191725163123, (11, 9) = -1.96363383725105, (11, 10) = -2.34025293598641, (11, 11) = -2.52909981686891, (11, 12) = -2.62024047869931, (11, 13) = -2.66301809488146, (11, 14) = -2.68268154837948, (11, 15) = -2.69157649927960, (11, 16) = -2.69554995185225, (11, 17) = -2.69730722992375, (11, 18) = -2.69807817758186, (11, 19) = -2.69841433230369, (11, 20) = -2.69856048527886, (11, 21) = -2.69862459845623, (11, 22) = -2.69865453167290, (11, 23) = -2.69867268287939, (11, 24) = -2.69869237852445, (11, 25) = -2.69872757167695, (11, 26) = -2.69880277010856, (11, 27) = -2.69897024989887, (11, 28) = -2.69934671724928, (11, 29) = -2.70019531250000, (12, 1) = 88.2002056515847, (12, 2) = 75.6725052605915, (12, 3) = 65.2163843281298, (12, 4) = 50.2957457752469, (12, 5) = 33.9778362938347, (12, 6) = 19.4889378828702, (12, 7) = 9.21145111328238, (12, 8) = 3.14668090532327, (12, 9) = 0.865015335380795e-2, (12, 10) = -1.48885415696429, (12, 11) = -2.17002830570858, (12, 12) = -2.47101994638412, (12, 13) = -2.60160334639512, (12, 14) = -2.65757048658506, (12, 15) = -2.68135582043060, (12, 16) = -2.69140318465282, (12, 17) = -2.69562852539876, (12, 18) = -2.69739963012287, (12, 19) = -2.69814032798427, (12, 20) = -2.69844990410614, (12, 21) = -2.69857998283484, (12, 22) = -2.69863653139994, (12, 23) = -2.69866541959327, (12, 24) = -2.69868944713954, (12, 25) = -2.69872638886708, (12, 26) = -2.69880229437275, (12, 27) = -2.69897006251682, (12, 28) = -2.69934665279911, (12, 29) = -2.70019531250000, (13, 1) = 212.942351062594, (13, 2) = 204.407863241860, (13, 3) = 192.852310449041, (13, 4) = 169.112345023827, (13, 5) = 129.709184621966, (13, 6) = 81.1275544412828, (13, 7) = 40.8611002373609, (13, 8) = 16.9812701000451, (13, 9) = 5.57967533249103, (13, 10) = .676047770238813, (13, 11) = -1.33996620211481, (13, 12) = -2.15392898726078, (13, 13) = -2.48042031376731, (13, 14) = -2.61116027872600, (13, 15) = -2.66353278518824, (13, 16) = -2.68453818355923, (13, 17) = -2.69297633500051, (13, 18) = -2.69637195950512, (13, 19) = -2.69774097819405, (13, 20) = -2.69829428533688, (13, 21) = -2.69851917836211, (13, 22) = -2.69861271212053, (13, 23) = -2.69865606573014, (13, 24) = -2.69868576547877, (13, 25) = -2.69872493730430, (13, 26) = -2.69880172276737, (13, 27) = -2.69896984150184, (13, 28) = -2.69934657781690, (13, 29) = -2.70019531250000, (14, 1) = 239.076755514012, (14, 2) = 236.081536444600, (14, 3) = 231.377517516101, (14, 4) = 220.073183560319, (14, 5) = 195.514343501259, (14, 6) = 149.755933550201, (14, 7) = 88.5760108265991, (14, 8) = 40.0127410081914, (14, 9) = 14.6812821286145, (14, 10) = 4.01218998170891, (14, 11) = -.143129764733769, (14, 12) = -1.72498492661880, (14, 13) = -2.32574215701891, (14, 14) = -2.55492675111780, (14, 15) = -2.64291819933079, (14, 16) = -2.67692134419163, (14, 17) = -2.69014146268037, (14, 18) = -2.69530980600786, (14, 19) = -2.69734057881215, (14, 20) = -2.69814249838062, (14, 21) = -2.69846133984409, (14, 22) = -2.69859056746230, (14, 23) = -2.69864754980966, (14, 24) = -2.69868247744579, (14, 25) = -2.69872366360231, (14, 26) = -2.69880122918730, (14, 27) = -2.69896965330897, (14, 28) = -2.69934651462986, (14, 29) = -2.70019531250000, (15, 1) = 242.787583178818, (15, 2) = 240.923505077373, (15, 3) = 237.803073889936, (15, 4) = 229.868632419040, (15, 5) = 211.146558598937, (15, 6) = 171.711129868456, (15, 7) = 109.599485934379, (15, 8) = 51.8362910328203, (15, 9) = 19.5660842529396, (15, 10) = 5.79120405458156, (15, 11) = .480594820388168, (15, 12) = -1.50742090421428, (15, 13) = -2.24937297699114, (15, 14) = -2.52785785730582, (15, 15) = -2.63322463639339, (15, 16) = -2.67341565955349, (15, 17) = -2.68886199142208, (15, 18) = -2.69483890414672, (15, 19) = -2.69716593067330, (15, 20) = -2.69807726797109, (15, 21) = -2.69843681897300, (15, 22) = -2.69858129502453, (15, 23) = -2.69864402434658, (15, 24) = -2.69868113037534, (15, 25) = -2.69872314674660, (15, 26) = -2.69880103063051, (15, 27) = -2.69896957817450, (15, 28) = -2.69934648954460, (15, 29) = -2.70019531250000, (16, 1) = 239.076755514012, (16, 2) = 236.081536444600, (16, 3) = 231.377517516101, (16, 4) = 220.073183560320, (16, 5) = 195.514343501259, (16, 6) = 149.755933550201, (16, 7) = 88.5760108265997, (16, 8) = 40.0127410081918, (16, 9) = 14.6812821286147, (16, 10) = 4.01218998170895, (16, 11) = -.143129764733754, (16, 12) = -1.72498492661880, (16, 13) = -2.32574215701891, (16, 14) = -2.55492675111780, (16, 15) = -2.64291819933079, (16, 16) = -2.67692134419163, (16, 17) = -2.69014146268037, (16, 18) = -2.69530980600786, (16, 19) = -2.69734057881215, (16, 20) = -2.69814249838062, (16, 21) = -2.69846133984409, (16, 22) = -2.69859056746230, (16, 23) = -2.69864754980966, (16, 24) = -2.69868247744579, (16, 25) = -2.69872366360231, (16, 26) = -2.69880122918730, (16, 27) = -2.69896965330897, (16, 28) = -2.69934651462986, (16, 29) = -2.70019531250000, (17, 1) = 212.942351062595, (17, 2) = 204.407863241860, (17, 3) = 192.852310449041, (17, 4) = 169.112345023828, (17, 5) = 129.709184621967, (17, 6) = 81.1275544412840, (17, 7) = 40.8611002373616, (17, 8) = 16.9812701000454, (17, 9) = 5.57967533249114, (17, 10) = .676047770238856, (17, 11) = -1.33996620211480, (17, 12) = -2.15392898726077, (17, 13) = -2.48042031376731, (17, 14) = -2.61116027872600, (17, 15) = -2.66353278518824, (17, 16) = -2.68453818355923, (17, 17) = -2.69297633500051, (17, 18) = -2.69637195950512, (17, 19) = -2.69774097819405, (17, 20) = -2.69829428533688, (17, 21) = -2.69851917836211, (17, 22) = -2.69861271212053, (17, 23) = -2.69865606573014, (17, 24) = -2.69868576547877, (17, 25) = -2.69872493730430, (17, 26) = -2.69880172276737, (17, 27) = -2.69896984150184, (17, 28) = -2.69934657781690, (17, 29) = -2.70019531250000, (18, 1) = 88.2002056515862, (18, 2) = 75.6725052605928, (18, 3) = 65.2163843281309, (18, 4) = 50.2957457752478, (18, 5) = 33.9778362938354, (18, 6) = 19.4889378828706, (18, 7) = 9.21145111328260, (18, 8) = 3.14668090532337, (18, 9) = 0.865015335385353e-2, (18, 10) = -1.48885415696427, (18, 11) = -2.17002830570857, (18, 12) = -2.47101994638411, (18, 13) = -2.60160334639511, (18, 14) = -2.65757048658506, (18, 15) = -2.68135582043060, (18, 16) = -2.69140318465282, (18, 17) = -2.69562852539876, (18, 18) = -2.69739963012287, (18, 19) = -2.69814032798427, (18, 20) = -2.69844990410614, (18, 21) = -2.69857998283484, (18, 22) = -2.69863653139994, (18, 23) = -2.69866541959327, (18, 24) = -2.69868944713954, (18, 25) = -2.69872638886708, (18, 26) = -2.69880229437275, (18, 27) = -2.69897006251682, (18, 28) = -2.69934665279911, (18, 29) = -2.70019531250000, (19, 1) = 15.2509102692802, (19, 2) = 12.6264291046441, (19, 3) = 10.5462806670589, (19, 4) = 7.66265417268544, (19, 5) = 4.61175757088344, (19, 6) = 1.94893912064135, (19, 7) = -0.655913599769431e-2, (19, 8) = -1.25191725163120, (19, 9) = -1.96363383725104, (19, 10) = -2.34025293598640, (19, 11) = -2.52909981686890, (19, 12) = -2.62024047869931, (19, 13) = -2.66301809488146, (19, 14) = -2.68268154837948, (19, 15) = -2.69157649927960, (19, 16) = -2.69554995185225, (19, 17) = -2.69730722992375, (19, 18) = -2.69807817758186, (19, 19) = -2.69841433230369, (19, 20) = -2.69856048527886, (19, 21) = -2.69862459845623, (19, 22) = -2.69865453167290, (19, 23) = -2.69867268287939, (19, 24) = -2.69869237852445, (19, 25) = -2.69872757167695, (19, 26) = -2.69880277010856, (19, 27) = -2.69897024989887, (19, 28) = -2.69934671724928, (19, 29) = -2.70019531250000, (20, 1) = .399782631544355, (20, 2) = -0.223230825286351e-1, (20, 3) = -.357156479975102, (20, 4) = -.825507840664786, (20, 5) = -1.32860914943847, (20, 6) = -1.78026685157638, (20, 7) = -2.12983346183127, (20, 8) = -2.36923077886600, (20, 9) = -2.51800244358535, (20, 10) = -2.60380139320853, (20, 11) = -2.65054483715776, (20, 12) = -2.67492220656519, (20, 13) = -2.68721033766193, (20, 14) = -2.69324013010359, (20, 15) = -2.69613566822986, (20, 16) = -2.69750183210721, (20, 17) = -2.69813712499983, (20, 18) = -2.69842905264433, (20, 19) = -2.69856199874022, (20, 20) = -2.69862239583270, (20, 21) = -2.69865047391492, (20, 22) = -2.69866531830100, (20, 23) = -2.69867716974714, (20, 24) = -2.69869424135636, (20, 25) = -2.69872834337068, (20, 26) = -2.69880308809385, (20, 27) = -2.69897037783243, (20, 28) = -2.69934676194892, (20, 29) = -2.70019531250000, (21, 1) = -2.16846353937697, (21, 2) = -2.23425128880022, (21, 3) = -2.28707572809903, (21, 4) = -2.36196107905337, (21, 5) = -2.44412091008287, (21, 6) = -2.52028432208616, (21, 7) = -2.58206991833073, (21, 8) = -2.62698442443030, (21, 9) = -2.65683978853650, (21, 10) = -2.67530598933600, (21, 11) = -2.68608729120433, (21, 12) = -2.69209644788596, (21, 13) = -2.69532209803656, (21, 14) = -2.69700104276086, (21, 15) = -2.69785291472636, (21, 16) = -2.69827602403542, (21, 17) = -2.69848245044229, (21, 18) = -2.69858169930828, (21, 19) = -2.69862895750863, (21, 20) = -2.69865157471866, (21, 21) = -2.69866311747724, (21, 22) = -2.69867077012410, (21, 23) = -2.69867951050463, (21, 24) = -2.69869524249446, (21, 25) = -2.69872876980332, (21, 26) = -2.69880326838836, (21, 27) = -2.69897045202630, (21, 28) = -2.69934678830801, (21, 29) = -2.70019531250000, (22, 1) = -2.60686981106158, (22, 2) = -2.61720825621594, (22, 3) = -2.62563727540733, (22, 4) = -2.63774342800024, (22, 5) = -2.65129960389903, (22, 6) = -2.66424133250227, (22, 7) = -2.67516577737618, (22, 8) = -2.68350330362916, (22, 9) = -2.68935662535048, (22, 10) = -2.69319142854833, (22, 11) = -2.69556414442440, (22, 12) = -2.69696419297891, (22, 13) = -2.69775825071771, (22, 14) = -2.69819388604856, (22, 15) = -2.69842625642532, (22, 16) = -2.69854727740723, (22, 17) = -2.69860905128715, (22, 18) = -2.69864009741297, (22, 19) = -2.69865562174525, (22, 20) = -2.69866364141685, (22, 21) = -2.69866853562029, (22, 22) = -2.69867318623206, (22, 23) = -2.69868058132789, (22, 24) = -2.69869571442844, (22, 25) = -2.69872897658248, (22, 26) = -2.69880335813610, (22, 27) = -2.69897048981909, (22, 28) = -2.69934680196416, (22, 29) = -2.70019531250000, (23, 1) = -2.68255270905730, (23, 2) = -2.68420269183551, (23, 3) = -2.68556927039621, (23, 4) = -2.68755569167303, (23, 5) = -2.68982174085213, (23, 6) = -2.69204188519929, (23, 7) = -2.69398007273551, (23, 8) = -2.69552020727773, (23, 9) = -2.69665146213758, (23, 10) = -2.69742910118661, (23, 11) = -2.69793453459984, (23, 12) = -2.69824780099437, (23, 13) = -2.69843426151571, (23, 14) = -2.69854146697746, (23, 15) = -2.69860129800890, (23, 16) = -2.69863385101707, (23, 17) = -2.69865119943290, (23, 18) = -2.69866033490219, (23, 19) = -2.69866522085018, (23, 20) = -2.69866814551416, (23, 21) = -2.69867062880846, (23, 22) = -2.69867415067918, (23, 23) = -2.69868102227875, (23, 24) = -2.69869591459725, (23, 25) = -2.69872906677682, (23, 26) = -2.69880339831608, (23, 27) = -2.69897050713020, (23, 28) = -2.69934680832502, (23, 29) = -2.70019531250000, (24, 1) = -2.69580300166118, (24, 2) = -2.69607068008258, (24, 3) = -2.69629574928807, (24, 4) = -2.69662647000063, (24, 5) = -2.69701006252744, (24, 6) = -2.69739450225433, (24, 7) = -2.69773989523952, (24, 8) = -2.69802385072690, (24, 9) = -2.69824051704782, (24, 10) = -2.69839565518402, (24, 11) = -2.69850083935603, (24, 12) = -2.69856887877845, (24, 13) = -2.69861113639229, (24, 14) = -2.69863647184190, (24, 15) = -2.69865120559918, (24, 16) = -2.69865955858620, (24, 17) = -2.69866421319792, (24, 18) = -2.69866682171521, (24, 19) = -2.69866840997105, (24, 20) = -2.69866969415176, (24, 21) = -2.69867137253520, (24, 22) = -2.69867450429699, (24, 23) = -2.69868118889116, (24, 24) = -2.69869599243597, (24, 25) = -2.69872910282171, (24, 26) = -2.69880341478773, (24, 27) = -2.69897051438702, (24, 28) = -2.69934681103519, (24, 29) = -2.70019531250000, (25, 1) = -2.69815364103587, (25, 2) = -2.69819775736826, (25, 3) = -2.69823537129639, (25, 4) = -2.69829118183365, (25, 5) = -2.69835687584068, (25, 6) = -2.69842403321791, (25, 7) = -2.69848588072185, (25, 8) = -2.69853822709065, (25, 9) = -2.69857949015361, (25, 10) = -2.69861008904377, (25, 11) = -2.69863160905534, (25, 12) = -2.69864606076830, (25, 13) = -2.69865538171212, (25, 14) = -2.69866118567008, (25, 15) = -2.69866469370480, (25, 16) = -2.69866676914122, (25, 17) = -2.69866799700861, (25, 18) = -2.69866877449387, (25, 19) = -2.69866940276336, (25, 20) = -2.69867019209140, (25, 21) = -2.69867161923136, (25, 22) = -2.69867462516171, (25, 23) = -2.69868124750408, (25, 24) = -2.69869602058751, (25, 25) = -2.69872911620663, (25, 26) = -2.69880342105693, (25, 27) = -2.69897051720914, (25, 28) = -2.69934681210573, (25, 29) = -2.70019531250000, (26, 1) = -2.69857559192165, (26, 2) = -2.69858297079499, (26, 3) = -2.69858934181473, (26, 4) = -2.69859887732956, (26, 5) = -2.69861024891079, (26, 6) = -2.69862207755826, (26, 7) = -2.69863320751553, (26, 8) = -2.69864286743489, (26, 9) = -2.69865069942731, (26, 10) = -2.69865668685300, (26, 11) = -2.69866103505642, (26, 12) = -2.69866405363369, (26, 13) = -2.69866606808988, (26, 14) = -2.69866736808164, (26, 15) = -2.69866818667951, (26, 16) = -2.69866870083484, (26, 17) = -2.69866904479461, (26, 18) = -2.69866933294639, (26, 19) = -2.69866969570361, (26, 20) = -2.69867034354260, (26, 21) = -2.69867169650278, (26, 22) = -2.69867466411125, (26, 23) = -2.69868126691862, (26, 24) = -2.69869603016244, (26, 25) = -2.69872912087587, (26, 26) = -2.69880342329623, (26, 27) = -2.69897051823817, (26, 28) = -2.69934681250191, (26, 29) = -2.70019531250000, (27, 1) = -2.69865211033863, (27, 2) = -2.69865336241077, (27, 3) = -2.69865445573850, (27, 4) = -2.69865610484201, (27, 5) = -2.69865809438942, (27, 6) = -2.69866019585749, (27, 7) = -2.69866221080965, (27, 8) = -2.69866399852749, (27, 9) = -2.69866548417810, (27, 10) = -2.69866665085402, (27, 11) = -2.69866752273294, (27, 12) = -2.69866814665308, (27, 13) = -2.69866857702353, (27, 14) = -2.69866886614376, (27, 15) = -2.69866906000736, (27, 16) = -2.69866919897384, (27, 17) = -2.69866932332963, (27, 18) = -2.69866948588351, (27, 19) = -2.69866977829260, (27, 20) = -2.69867038746753, (27, 21) = -2.69867171953946, (27, 22) = -2.69867467603810, (27, 23) = -2.69868127301993, (27, 24) = -2.69869603324795, (27, 25) = -2.69872912241713, (27, 26) = -2.69880342405215, (27, 27) = -2.69897051859236, (27, 28) = -2.69934681264019, (27, 29) = -2.70019531250000, (28, 1) = -2.69866603047084, (28, 2) = -2.69866625169965, (28, 3) = -2.69866644691331, (28, 4) = -2.69866674348742, (28, 5) = -2.69866710514694, (28, 6) = -2.69866749259564, (28, 7) = -2.69866787060424, (28, 8) = -2.69866821286631, (28, 9) = -2.69866850389706, (28, 10) = -2.69866873830742, (28, 11) = -2.69866891844754, (28, 12) = -2.69866905155222, (28, 13) = -2.69866914728661, (28, 14) = -2.69866921625670, (28, 15) = -2.69866926986058, (28, 16) = -2.69866932201265, (28, 17) = -2.69866939402006, (28, 18) = -2.69866952574703, (28, 19) = -2.69866980038930, (28, 20) = -2.69867039952334, (28, 21) = -2.69867172602142, (28, 22) = -2.69867467947623, (28, 23) = -2.69868127482055, (28, 24) = -2.69869603417946, (28, 25) = -2.69872912289265, (28, 26) = -2.69880342429014, (28, 27) = -2.69897051870583, (28, 28) = -2.69934681268504, (28, 29) = -2.70019531250000, (29, 1) = -2.69866813533281, (29, 2) = -2.69866820985674, (29, 3) = -2.69866827617016, (29, 4) = -2.69866837749853, (29, 5) = -2.69866850212862, (29, 6) = -2.69866863715926, (29, 7) = -2.69866877073077, (29, 8) = -2.69866889363839, (29, 9) = -2.69866900008058, (29, 10) = -2.69866908759731, (29, 11) = -2.69866915646738, (29, 12) = -2.69866920891275, (29, 13) = -2.69866924843909, (29, 14) = -2.69866927960774, (29, 15) = -2.69866930859281, (29, 16) = -2.69866934517168, (29, 17) = -2.69866940758526, (29, 18) = -2.69866953354284, (29, 19) = -2.69866980479120, (29, 20) = -2.69867040196861, (29, 21) = -2.69867172735934, (29, 22) = -2.69867468019801, (29, 23) = -2.69868127520480, (29, 24) = -2.69869603438140, (29, 25) = -2.69872912299730, (29, 26) = -2.69880342434324, (29, 27) = -2.69897051873145, (29, 28) = -2.69934681269525, (29, 29) = -2.70019531250000}, datatype = float[8], storage = rectangular, order = Fortran_order, shape = []):

``

Convert to a "spline"

 

``

SolSpline := GridToSpline(SampleGrid, TCut, RCut, nt, nr): 

``

plot3d(SolSpline, -TCut .. TCut, 0 .. RCut)

 

 

My Integrals are not straight forward. I have some analytic function

 

 

Vapprox := proc (rho) options operator, arrow; (-1)*1.6583067*10^8*((1/175)*rho-5/7)^10+(-1)*4.8124279*10^8*((1/175)*rho-5/7)^11+1.8264850*10^8*((1/175)*rho-5/7)^12+(-1)*1.0974043*10^7*((1/175)*rho-5/7)^5+(-1)*3.2006497*10^7*((1/175)*rho-5/7)^6+(-1)*1.3012638*10^7*((1/175)*rho-5/7)^7+9.8671189*10^7*((1/175)*rho-5/7)^8+1.8035498*10^8*((1/175)*rho-5/7)^9+3.4504836*10^7+(-1)*8.0188022*10^7*((1/175)*rho-5/7)^2+4.3558999*10^7*((1/175)*rho-5/7)^3+8.4893913*10^7*((1/175)*rho-5/7)^4+5.6273392*10^7*((1/175)*rho-5/7)^16+8.8592685*10^7*((1/175)*rho-5/7)^17+(-1)*1.0795555*10^7*((1/175)*rho-5/7)^18+5.9763417*10^8*((1/175)*rho-5/7)^13+(-1)*1.3153369*10^8*((1/175)*rho-5/7)^14+(-1)*3.6494514*10^8*((1/175)*rho-5/7)^15+(-1)*312799.54*rho end proc

Vapprox := proc (rho) options operator, arrow; (-1)*1.6583067*100000000*(((1/175)*rho-5/7)^10)+(-1)*4.8124279*100000000*(((1/175)*rho-5/7)^11)+1.8264850*100000000*(((1/175)*rho-5/7)^12)+(-1)*1.0974043*10000000*(((1/175)*rho-5/7)^5)+(-1)*3.2006497*10000000*(((1/175)*rho-5/7)^6)+(-1)*1.3012638*10000000*(((1/175)*rho-5/7)^7)+9.8671189*10000000*(((1/175)*rho-5/7)^8)+1.8035498*100000000*(((1/175)*rho-5/7)^9)+3.4504836*10000000+(-1)*8.0188022*10000000*(((1/175)*rho-5/7)^2)+4.3558999*10000000*(((1/175)*rho-5/7)^3)+8.4893913*10000000*(((1/175)*rho-5/7)^4)+5.6273392*10000000*(((1/175)*rho-5/7)^16)+8.8592685*10000000*(((1/175)*rho-5/7)^17)+(-1)*1.0795555*10000000*(((1/175)*rho-5/7)^18)+5.9763417*100000000*(((1/175)*rho-5/7)^13)+(-1)*1.3153369*100000000*(((1/175)*rho-5/7)^14)+(-1)*3.6494514*100000000*(((1/175)*rho-5/7)^15)+(-1)*312799.54*rho end proc

 

````

Vapprox(2);

-208375.744

 

 

And I need to take numeric derivatives in the t and r directions. I also need to make sure only numeric arguements are passed, so I've built the following. (Which have global variables, not elegant I know).

 

SSpline := proc (t, r) if not (t::numeric and r::numeric) then return ('procname')(args) end if; SolSpline(t, r) end proc:
 

 

``

dtSpline := proc (tval, r) if not (tval::numeric and r::numeric) then return ('procname')(args) end if; return fdiff(SSpline(t, r), t = tval) end proc:

NULL

drSpline := proc (t, rval) if not (t::numeric and rval::numeric) then return ('procname')(args) end if; return fdiff(SSpline(t, r), r = rval) end proc:

``

````

Now my Integrals.

 

These only take 50 seconds when I am using a function outputted for dsolve as my "SSpline" function, however with the ArrayInterpolation it takes more time as you can see...

 

 

rLength := RCut;

rLength := .3

 

 

 

with(CodeTools):

``

A1 := Usage(evalf(4*Pi*(int(r^2*Vapprox(SSpline(t, r)), [r = 0 .. RCut, t = -rLength .. rLength], digits = 7, method = _CubaCuhre))))

memory used=97.96GiB, alloc change=0 bytes, cpu time=5.59m, real time=4.80m, gc time=114.42s
A1 := -20234.83492

 

``

A3 := CodeTools:-Usage(evalf(4*Pi*(int(.5*r^2*dtSpline(t, r)^2, [r = 0 .. rLength, t = -rLength .. rLength], digits = 7, method = _CubaCuhre))))

Warning,  computation interrupted

 

A4 := CodeTools:-Usage(evalf(4*Pi*(int(.5*r^2*drSpline(t, r)^2, [r = 0 .. rLength, t = -rLength .. rLength], digits = 7, method = _CubaCuhre))))

Warning,  computation interrupted

 

A1+A3+A4

 

``

``

``


Download SimplifiedIntegration.mwSimplifiedIntegration.mw

@Carl Love 

Thanks for the quick reply. I would like to test this out, however I'm using Maple 18.02, which doesn't seem to incorporate the lowess function. :(

@ecterrab 

Perfect. 

I hadn't tried redefining the g_[] metric. In hindsight it seems an obvious solution. 

This is the solution I was chasing.

Thanks

@acer  

I've had a look at this problem more closely. The issues lies with my stucture of V. In the rho range at hand, I get complex exponentials which are not dealt with correctly. Fundamentally this is an issue with the V function I have defined. I will need to correct this and/or work out how to handle what I am expecting in this complex regime. 

 

All the numerical solving techniques you have provided will then be very useful. Thanks again for your help and hopefully this thread helps others in the future. 

 

Cheers

@acer Yes, I was busy trying to get the numerical ODE to solve and hadn't payed attention to the actual IVP I had sent. I'm not sure if you're interested in the problem, but its a physics one. 

I have the function V, which is a potential. It contains three parts V_0, VCWc and V_TDD. Because of the complications of V_TDD, I have given it default values of lambda and mu (defined at the start of the sheet). The other two, I have used more accurate coefficients for lambda and mu, namely the solutions of the fsolve argument where I've set rho=250 and 125^2. Func is essentially V_C + VCW (if mu and lambda were undefined variables). Lets call this "Func". I take the first derivative of of Func wrt to rho and at the point rho=246 this derivative should be zero. Similarly the second derivative should be 125^2. These conditions give me values of mu and lambda as a function of kappa. 

The physical problem is that I am writting a modified potential to the Higgs field. Its derivative at the known minimum (electroweak minumum 246 Gev) should be zero and the second derivative should be the higgs mass, 125^2. 

In regards to accuracy, I'm trying to solve this with what ever accuracy I can get. The ODE I'm solving is like a ball rolling down a hill with friction. The potential looks like 

where the blue line is the negative potential. I need to find the initial rho value which gives a result such that it rolls down the potential and lands onto the second bump. I do this via a bisection. Essentially I pick some initial point p(0.001)=190, solve my ODE see if it goes over or under the bump. Then repeat with a higher/lower value until I narrow down to some degree of accuracy where it sits on the second hill for a little bit before falling back to its minimum (in this case around 0). I'm using the Student[NumericalAnalysis]-Roots procedure to do this. 

This means, I need to repetitively solve this IVP and require speed. I spose I dont need a rediculous level of accuracy for each iteration. 

I will have a play with some polynomial fitting functions to approximate this derivative. 

Thanks again!

 

 

@acer 

Hi Acer, thank you VERY much for this amazing response. I've learnt a lot about useful methods for numerical differentiion. 

I just have one further, hopefully small problem. 

I had intentionally used the negative derviative (was used for an entirely different reason), however in this problem, it needs to be positive. The solution I'm looking for ranges over a larger range of rho. Typically for Rho from -50 to 250. 

The solutions presented here work very well for rho from 1 to 10. However the symbolic differentiation breaks down between rho = 50 and 150. (this was one of the main reasons I chose to numerically differentiate). 

The solutions you have presented for mapping a polynomial of degree n is very useful, however I am having difficulty mapping a polynomial (ideally to the dVdrho function you have created, or alternatively to the fdiff function). However I am unable to do this, the mapping to the fdiff function seems to take an unpractically long time, do you know of a method to speed this mapping up?. 

Thanks in advance. 

I've attached your previous worksheet with a few comments on the end. 

 


``

``

restart

plots:-setoptions(gridlines = false);

Eq := diff(rho(r), r, r)+2*(diff(rho(r), r))/r-dV/(d*rho(r));

diff(diff(rho(r), r), r)+2*(diff(rho(r), r))/r-dV/(d*rho(r))

(1)

m[h] := 125:

(2)

(3)

y[t] := sqrt(2)*m[t]/v:

kappabare := -m[h]^2*kappa/v: 

mu := sqrt(1/2*(v*kappabare+m[h]^2)):``

lambda := (-v*kappabare+m[h]^2)/(2*v^2):

n := [4, 2, 2, 1, -12, 1, 1]:NULL

mz := ((g[1]^2+g[2]^2)*(1/4))*rho^2:

delta := (1/4)*(mz+(11/6*(g[1]^2+g[2]^2))*T^2)^2-(11*g[1]^2*g[2]^2*(1/12))*T^2*((11/3)*T^2+rho^2):

ms := [(1/4)*g[2]^2*rho^2, (1/4)*g[2]^2*rho^2+(11/6)*g[2]^2*T^2, mz, (1/2)*mz+(11/12*(g[1]^2+g[2]^2))*T^2+(1/2)*sqrt(delta), (1/2)*y[t]^2*rho^2, 3*lambda*rho^2+2*kappabare*rho-mu^2+(1/8)*g[2]^2*T^2+(1/16*(g[1]^2+g[2]^2))*T^2+(1/4)*y[t]^2*T^2, (1/2)*mz+(11/12*(g[1]^2+g[2]^2))*T^2-(1/2)*sqrt(delta)]:

s := [1, 1, 1, 1, 1/2, 0, 1]:

Vtemp := 0;

0

(4)

for i to 7 do Vtemp := Vtemp+T^4*n[i]*(Int(Re(x^2*ln(1+(-1)^(2*s[i]+1)*exp(-sqrt(x^2+ms[i]/T^2)))), x = 0 .. infinity))/(2*Pi^2) end do:


The integrals can be combined into a single (semi-infinite) integral. I didn't look to see whether the resulting single integrand could be better simplified.

Vtemp := simplify(IntegrationTools:-Combine(Vtemp)):

V_TDD := unapply(Vtemp2, rho, T, kappa);

proc (rho, T, kappa) options operator, arrow; Int(0.5066059180e-1*T^4*x^2*(ln(abs(-1.+exp(-1.*((x^2*T^2+.5016825145*T^2+0.6841125200e-1*rho^2+.5000000000*(0.4680099400e-2*rho^4+0.2043678598e-1*rho^2*T^2+0.7493488170e-1*T^4)^(1/2))/T^2)^(1/2))))+ln(abs(-1.+exp(-1.*(-(-x^2*T^2-.5016825145*T^2-0.6841125200e-1*rho^2+.5000000000*(0.4680099400e-2*rho^4+0.2043678598e-1*rho^2*T^2+0.7493488170e-1*T^4)^(1/2))/T^2)^(1/2))))+ln(abs(-1.+exp(-1.*((x^2*T^2+.3872942693*rho^2*kappa+.3872942693*rho^2-127.0325203*kappa*rho+7812.500000*kappa-7812.500000+.3353301441*T^2)/T^2)^(1/2))))+2.*ln(abs(-1.+exp(-1.*((x^2*T^2+.1057397562*rho^2+.7754248784*T^2)/T^2)^(1/2))))+4.*ln(abs(-1.+exp(-1.*((x^2*T^2+.1057397562*rho^2)/T^2)^(1/2))))-12.*ln(abs(1.+exp(-1.*((x^2*T^2+.4965092801*rho^2)/T^2)^(1/2))))+2.*ln(abs(-1.+exp(-1.*((x^2*T^2+.1368225040*rho^2)/T^2)^(1/2))))), x = 0. .. Float(infinity), digits = 15, epsilon = 0.1e-11, method = _d01amc) end proc

(5)


We can symbolically differentiate the expression Vtemp wrt rho, before unapplying it. This gets used in the latter two of the four approaches below to solve the IVP.

dVtemp := simplify(combine(diff(Vtemp, rho)), size):

mu := 'mu':

func := -Typesetting:-delayDotProduct(1/2, m, true)*rho^2-(15625/738)*kappa*rho^3+(1/4)*lambda*rho^4+0.1048209004e-2*rho^4*(ln(0.1747302469e-5*rho^2-1/10000000000)-3/2)/Pi^2+0.8775186375e-3*rho^4*(ln(0.2260931060e-5*rho^2-1/10000000000)-3/2)/Pi^2-0.4622277472e-1*rho^4*(ln(0.8204595150e-5*rho^2-1/10000000000)-3/2)/Pi^2+(1/64)*(3*lambda*rho^2-(15625/123)*kappa*rho-m)^2*(ln((1/20172)*lambda*rho^2-(15625/7443468)*kappa*rho-Typesetting:-delayDotProduct(1/60516, m, true)-1/10000000000)-3/2)/Pi^2:
NULL

eq1 := subs(rho = 246, diff(func, rho)) = 0:

eq1p := unapply(eq1, kappa):

eq2 := subs(rho = 246, diff(func, rho, rho)) = 125^2:``

eq2p := unapply(eq2, kappa):


Make cmu and clambda return unevaluated for nonnumeric arguments, so that we can hit VCWc_alt with the D command.

Give them option remember,system so that their results are memoized (up until garbage collection). This might help speed, since VCWc_alt contains two identical calls to each of these procedures. A further refinement would be to notice that the fsolve call is the same, four times, and there should be no need to have each of clu and clambda call fsolve with the same arguments just to be able to pick off two seperate variables' results. (I didn't go that far...)

It's possible that `func` (or similar) could be used as an expression, simplifed to hopefully reduce roundoff error in the float evaluation of the complicated expression, and then used in VCWc. (I didn't do that...)

cmu := proc (kappa) options remember, system; if not type(kappa, numeric) then return ('procname')(args) end if; rhs(fsolve({eq1p(kappa), eq2p(kappa)})[2]) end proc:

cmu(kappa)

(6)

(7)

clambda := proc (kappa) options remember, system; if not type(kappa, numeric) then return ('procname')(args) end if; rhs(fsolve({eq1p(kappa), eq2p(kappa)})[1]) end proc:

clambda(kappa)

(8)

Now to add the last two parts of the V function. (Here I use the actual values of cmu and clambda directly, it helps in some circumstances)

NULL

VCWc := proc (rho, kappa) options operator, arrow; 0.1048209004e-2*rho^4*(ln(0.1747302469e-5*rho^2-1/10000000000)-3/2)/Pi^2+0.8775186375e-3*rho^4*(ln(0.2260931060e-5*rho^2-1/10000000000)-3/2)/Pi^2+(-1)*0.4622277472e-1*rho^4*(ln(0.8204595150e-5*rho^2-1/10000000000)-3/2)/Pi^2+(1/64)*(3*rhs(fsolve({eq1p(kappa), eq2p(kappa)})[1])*rho^2-(15625/123)*kappa*rho-rhs(fsolve({eq1p(kappa), eq2p(kappa)})[2]))^2*(ln((1/20172)*rhs(fsolve({eq1p(kappa), eq2p(kappa)})[1])*rho^2-(15625/7443468)*kappa*rho-(1/60516)*rhs(fsolve({eq1p(kappa), eq2p(kappa)})[2])-1/10000000000)-3/2)/Pi^2 end proc:
 


Now that clambda and cmu return unevaluated for nonnumeric argument we can successfully
differentiate the procedure altVCWc wrt its first argument (rho), using the D command.

altVCWc := proc (rho, kappa) options operator, arrow; 0.1048209004e-2*rho^4*(ln(0.1747302469e-5*rho^2-1/10000000000)-3/2)/Pi^2+0.8775186375e-3*rho^4*(ln(0.2260931060e-5*rho^2-1/10000000000)-3/2)/Pi^2+(-1)*0.4622277472e-1*rho^4*(ln(0.8204595150e-5*rho^2-1/10000000000)-3/2)/Pi^2+(1/64)*(3*clambda(kappa)*rho^2-(15625/123)*kappa*rho-cmu(kappa))^2*(ln((1/20172)*clambda(kappa)*rho^2-(15625/7443468)*kappa*rho-(1/60516)*cmu(kappa)-1/10000000000)-3/2)/Pi^2 end proc:


We expect these two to give very similar results.

plot('Re(VCWc(rho, 1.5))'-'Re(altVCWc(rho, 1.5))', rho = .1 .. 20, view = -0.1e-5 .. 0.1e-5, size = [600, 100], thickness = 3, adaptive = false, numpoints = 20);

 

dVCWc := D[1](altVCWc);

proc (rho, kappa) options operator, arrow; 0.4192836016e-2*rho^3*(ln(0.1747302469e-5*rho^2-1/10000000000)-3/2)/Pi^2+0.3663076361e-8*rho^5/((0.1747302469e-5*rho^2-1/10000000000)*Pi^2)+0.3510074550e-2*rho^3*(ln(0.2260931060e-5*rho^2-1/10000000000)-3/2)/Pi^2+0.3968018287e-8*rho^5/((0.2260931060e-5*rho^2-1/10000000000)*Pi^2)-.1848910989*rho^3*(ln(0.8204595150e-5*rho^2-1/10000000000)-3/2)/Pi^2-0.7584783066e-6*rho^5/((0.8204595150e-5*rho^2-1/10000000000)*Pi^2)+(1/32)*(3*clambda(kappa)*rho^2-(15625/123)*kappa*rho-cmu(kappa))*(ln((1/20172)*clambda(kappa)*rho^2-(15625/7443468)*kappa*rho-(1/60516)*cmu(kappa)-1/10000000000)-3/2)*(6*clambda(kappa)*rho-(15625/123)*kappa)/Pi^2+(1/64)*(3*clambda(kappa)*rho^2-(15625/123)*kappa*rho-cmu(kappa))^2*((1/10086)*clambda(kappa)*rho-(15625/7443468)*kappa)/(((1/20172)*clambda(kappa)*rho^2-(15625/7443468)*kappa*rho-(1/60516)*cmu(kappa)-1/10000000000)*Pi^2) end proc

(9)

V_0c := proc (rho, kappa) options operator, arrow; -(1/2)*cmu(kappa)*rho^2-(1/3)*m[h]^2*kappa*rho^3/v+(1/4)*clambda(kappa)*rho^4 end proc;

proc (rho, kappa) options operator, arrow; -(1/2)*cmu(kappa)*rho^2-(1/3)*m[h]^2*kappa*rho^3/v+(1/4)*clambda(kappa)*rho^4 end proc

(10)

dV_0c := D[1](V_0c);

proc (rho, kappa) options operator, arrow; -cmu(kappa)*rho-m[h]^2*kappa*rho^2/v+clambda(kappa)*rho^3 end proc

(11)

V := proc (rho, T, kappa) Digits := 15; V_0c(rho, kappa)+Re(VCWc(rho, kappa))+V_TDD(rho, T, kappa) end proc;

proc (rho, T, kappa) Digits := 15; V_0c(rho, kappa)+Re(VCWc(rho, kappa))+V_TDD(rho, T, kappa) end proc

(12)

 


Let's also build dVdrho which instead uses the "derivatives" of V_0c, VCWc, and V_TDD. This provides ODEProc_alt below which does not need fdiff or any other numeric differentiation scheme.

dVdrho := subs([V_0c = dV_0c, VCWc = dVCWc, V_TDD = dV_TDD], proc (rho, T, kappa) Digits := 15; V_0c(rho, kappa)+Re(VCWc(rho, kappa))+V_TDD(rho, T, kappa) end proc);

proc (rho, T, kappa) Digits := 15; dV_0c(rho, kappa)+Re(dVCWc(rho, kappa))+dV_TDD(rho, T, kappa) end proc

(13)

V(20, 100, 1.5);

301107.785317594+2667.1610609935/Pi^2+Int(5066059.18000000*x^2*(ln(abs(-1.+exp(-1.*(x^2+.642041492984909)^(1/2))))+ln(abs(-1.+exp(-1.*(x^2+.366796436175091)^(1/2))))+ln(abs(-1.+exp(-1.*(x^2+.383587010130000)^(1/2))))+2.*ln(abs(-1.+exp(-1.*(x^2+.779654468648000)^(1/2))))+4.*ln(abs(-1.+exp(-1.*(x^2+0.422959024800000e-2)^(1/2))))-12.*ln(abs(1.+exp(-1.*(x^2+0.198603712040000e-1)^(1/2))))+2.*ln(abs(-1.+exp(-1.*(x^2+0.547290016000000e-2)^(1/2))))), x = 0. .. Float(infinity), digits = 15, epsilon = 0.1e-11, method = _d01amc)

(14)

evalf(V(20, 100, 1.5));

-227049562.0

(15)

CodeTools:-Usage(plot(V(rho, 100, 1.5), rho = 0 .. 20, adaptive = false, numpoints = 30, size = [600, 100]))

 


I suppose that you realize that by using a difference formula like (f(x-h) - f(x))/h you are obtaining the negative of the derivative. If you flip the sign back then your DE might encounter a singularity early on. I presume you have it the way you intend it.

VDiff := proc (rho, T, kappa) options operator, arrow; evalf(V(rho+0.1e-1, T, kappa)-V(rho, T, kappa))/0.1e-1 end proc;

proc (rho, T, kappa) options operator, arrow; evalf(V(rho+0.1e-1, T, kappa)-V(rho, T, kappa))/0.1e-1 end proc

(16)

 

ODEProc := proc (rho) local t; if not rho::numeric then return ('procname')(args) end if; Digits := 12; fdiff(V(t, 100, 1.5), t = rho) end proc:

ODEProc_alt := proc (rho) if not type(rho, numeric) then return ('procname')(args) end if; Digits := 12; evalf(dVdrho(rho, 100, 1.5)) end proc:

 

 

Here its clear the symbolic differentiation breaks down, however the fdiff function seems to work reasonably.

 

plot([ODEProc_alt(rho), ODEProc(rho)], rho = 0 .. 250, style = [line, point], adaptive = false, numpoints = 30, size = [600, 100]);

 

 

 

Digits := 12; CodeTools:-Usage(numapprox:-chebpade(ODEProc_alt(rho), rho = 1 .. 250, [6, 4])); vdiffapprox_alt := eval(%, T = orthopoly[T]); Digits := 10

12

(17)

``The same code applied to the fdiff function seems to take an unpractically long time.

 

 

``


Download help2.mw

 

 

@acer 

Hi acer, thanks for the indepth response. For some reason I cant enter the contents of this maple file, but here is a link. I've added more complexity (probably too much), but my problem still seems trivial, in that I need to get dsolve to run over some black box function. 

Help.mw

 

Thanks!

@Kitonum @tomleslie 

 

A simplified version of the function is essentially this guy:

 

g:=(z)->Int(Re(x^2*ln(1+exp(-sqrt(x^2+(-.23*z^2)*(1/45000))))), x = 0 .. infinity, digits = 10, method = _d01amc)

I'm using this specific integration technique as its the fastest for solving the full form of the equation. 

@Kitonum, I've tried using fsolve. I can plot the function, and see roughly where the roots lie. Running fsolve(deriv) finds roots which are not accurate and in some cases are blatently wrong.

However you've solved my problem. It was using the ' ' - to wrap around my function. I can now use the Roots package and it gives the correct solutions. Perfect. 

Is there page you could link me to, to help me understand the exact meaning of the ' '?

 

Thanks guys!!

 

 

 

 

 

@ecterrab 

Perhaps I am misunderstanding the math. I am trying to define some tensor e[mu,~a] as you have done above in your second line. I would expect the components of e[mu,~a] to be the same as defined? (Not contracted with the metric). 

The scenario I am picturing is illustrated in the following code (I'm using the same example as before as I've already got the maple code :)). 


restart; with(Physics); with(Tetrads); Coordinates(X = [t, r, theta, z])

LE := -dr^2-Physics:-`*`(r^2, dtheta^2)-dz^2+dt^2;

Setup(mathematicalnotation = true, metric = LE):

Setup(tetradmetric = Matrix(4, `<,>`(1, -1, -1, -1), shape = diagonal));

e[a, `~mu`] = Matrix(4, `<,>`(1, 1, 1/r, 1), shape = diagonal)

 

Lets say I define some Tetrad as below. 

 

e[a, `~mu`] = (Matrix(4, 4, {(1, 1) = 1, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (2, 1) = 0, (2, 2) = 1, (2, 3) = 0, (2, 4) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 1/r, (3, 4) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 1}))

(4)

Define(%, quiet):

e[a, `~mu`, matrix]

 

Its components are what I expect. 

 

e[a, `~mu`] = (Matrix(4, 4, {(1, 1) = 1, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (2, 1) = 0, (2, 2) = 1, (2, 3) = 0, (2, 4) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 1/r, (3, 4) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 1}))

(5)

 

The invserse Tetrad is also what I would expect. 

 

e[`~a`, mu, matrix]

e[`~a`, mu] = (Matrix(4, 4, {(1, 1) = 1, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (2, 1) = 0, (2, 2) = 1, (2, 3) = 0, (2, 4) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = r, (3, 4) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 1}))

(6)

Define(clear, e):

`Defined as tensors`

(7)

 

What if, I only knew the inverse Tetrad. I would like to be able to input the Inverse Tetrad without having to manually convert it to the "non-inverted" form. I would expect to define it as follows.

``

e[`~a`, mu] = (Matrix(4, 4, {(1, 1) = 1, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (2, 1) = 0, (2, 2) = 1, (2, 3) = 0, (2, 4) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = r, (3, 4) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 1}))

e[`~a`, mu] = (Matrix(4, 4, {(1, 1) = 1, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (2, 1) = 0, (2, 2) = 1, (2, 3) = 0, (2, 4) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = r, (3, 4) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 1}))

(8)

Define(%);

{Physics:-D_[mu], Physics:-Dgamma[mu], Physics:-Psigma[mu], Physics:-Ricci[mu, nu], Physics:-Riemann[mu, nu, alpha, beta], Physics:-Weyl[mu, nu, alpha, beta], Physics:-d_[mu], e[`~a`, mu], Physics:-Tetrads:-e_[a, mu], Physics:-Tetrads:-eta_[a, b], Physics:-g_[mu, nu], Physics:-Tetrads:-gamma_[a, b, c], Physics:-Tetrads:-l_[mu], Physics:-Tetrads:-lambda_[a, b, c], Physics:-Tetrads:-m_[mu], Physics:-Tetrads:-mb_[mu], Physics:-Tetrads:-n_[mu], Physics:-Christoffel[mu, nu, alpha], Physics:-Einstein[mu, nu], Physics:-KroneckerDelta[mu, nu], Physics:-LeviCivita[alpha, beta, mu, nu], Physics:-SpaceTimeVector[mu](X)}

(9)

e[`~a`, mu, matrix];

e[`~a`, mu] = (Matrix(4, 4, {(1, 1) = 1, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (2, 1) = 0, (2, 2) = 1, (2, 3) = 0, (2, 4) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = r^3, (3, 4) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 1}))

(10)

e[a, `~mu`, matrix]

e[a, `~mu`] = (Matrix(4, 4, {(1, 1) = 1, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (2, 1) = 0, (2, 2) = 1, (2, 3) = 0, (2, 4) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = r, (3, 4) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 1}))

(11)

``Its obviously now not the Tetrad I started with. I have a feeling I'm being naive here somewhere, perhaps you could point out my error in logic/maple. 


Download InverseTetradTensors.mw

@ecterrab 

Hi Edgardo, 

Thanks for the tips. I'm not sure if you want me to keep posting bugs, or just wait a couple of months for the package to be ironed out. 

I'll add just one last one. 

Defining Tensors with tetrad indicies - It doesn't like the ordering or contravariant definitions. 

Eg.

e[mu,~a] = <some matrix>; Define(%)

e[mu,~a,matrix] != <some matrix> but

e[a,~mu, matrix] = <some matrix>

Also

e[~a,mu] = <some matrix>; Define(%)

e[~a,mu,matrix] != <some matrix>

e[~mu,a,matrix] = <some matrix>

In these example's i'm using a non-minkowskian metric. 

@ecterrab 

Hi, I've noticed some odd calculations of Tensors with tetrad indices. Here is a simple example. 



restart

with(Physics):

with(Physics:-Tetrads)

[e_, eta_, gamma_, l_, lambda_, m_, mb_, n_]

(1)

``

Use Cylindrical Co-ords as an example.

 

LE := -dr^2-Physics:-`*`(r^2, dtheta^2)-dz^2+dt^2;

-dtheta^2*r^2-dr^2+dt^2-dz^2

(2)

``

[t, r, theta, z]

(3)

Coordinates(X = [t, r, theta, z]):

`Systems of spacetime Coordinates are: `*{X = (t, r, theta, z)}

(4)

Setup(mathematicalnotation = true, metric = LE):

Setup(tetradmetric = Matrix(4, `<,>`(1, -1, -1, -1), shape = diagonal));

[tetradmetric = {(1, 1) = 1, (2, 2) = -1, (3, 3) = -1, (4, 4) = -1}]

(5)

``

 

SumOverRepeatedIndices(Physics:-`*`(Physics:-`*`(g_[`~mu`, `~nu`], e_[`~a`, mu]), e_[`~b`, nu]))

Physics:-Tetrads:-e_[`~a`, 1]*Physics:-Tetrads:-e_[`~b`, 1]-Physics:-Tetrads:-e_[`~a`, 2]*Physics:-Tetrads:-e_[`~b`, 2]-Physics:-Tetrads:-e_[`~a`, 3]*Physics:-Tetrads:-e_[`~b`, 3]/r^2-Physics:-Tetrads:-e_[`~a`, 4]*Physics:-Tetrads:-e_[`~b`, 4]

(6)

TensorArray(%)

Matrix(4, 4, {(1, 1) = 1, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (2, 1) = 0, (2, 2) = -1, (2, 3) = 0, (2, 4) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = -1, (3, 4) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = -1})

(7)

 

Above works as expected. What about a custom Tetrad?

e[a, `~mu`] = Matrix(4, `<,>`(1, 1, 1/r, 1), shape = diagonal)

e[a, `~mu`] = (Matrix(4, 4, {(1, 1) = 1, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (2, 1) = 0, (2, 2) = 1, (2, 3) = 0, (2, 4) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 1/r, (3, 4) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 1}))

(8)

``

Define(%, quiet):

e[`~a`, mu, matrix]

e[`~a`, mu] = (Matrix(4, 4, {(1, 1) = 1, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (2, 1) = 0, (2, 2) = 1, (2, 3) = 0, (2, 4) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = r, (3, 4) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 1}))

(9)

 

Inverse is as expected. Lets perform the same summation as above. 

 

 

SumOverRepeatedIndices(Physics:-`*`(Physics:-`*`(g_[`~mu`, `~nu`], e[`~a`, mu]), e[`~b`, nu]))

e[1, `~a`]*e[1, `~b`]-e[2, `~a`]*e[2, `~b`]-e[3, `~a`]*e[3, `~b`]/r^2-e[4, `~a`]*e[4, `~b`]

(10)

TensorArray(%);

Matrix(4, 4, {(1, 1) = 1, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (2, 1) = 0, (2, 2) = -1, (2, 3) = 0, (2, 4) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = -1/r^4, (3, 4) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = -1})

(11)

%[3, 3]

-1/r^4

(12)

 

This guy appears? I assume its not handling tetrad indices correctly. 

 

 

T2[`~a`, `~b`] = e[1, `~a`]*e[1, `~b`]-e[2, `~a`]*e[2, `~b`]-e[3, `~a`]*e[3, `~b`]/r^2-e[4, `~a`]*e[4, `~b`]

T2[`~a`, `~b`] = e[1, `~a`]*e[1, `~b`]-e[2, `~a`]*e[2, `~b`]-e[3, `~a`]*e[3, `~b`]/r^2-e[4, `~a`]*e[4, `~b`]

(13)

Define(%, quiet):

T2[`~a`, `~b`, matrix]

T2[`~a`, `~b`] = (Matrix(4, 4, {(1, 1) = 1, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (2, 1) = 0, (2, 2) = -1, (2, 3) = 0, (2, 4) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = -1, (3, 4) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = -1}))

(14)

 

In Matrix form, when defined as a tensor, it doesn't appear.

 

 

T2[`~3`, `~3`]

-1/r^4

(15)

 

But appears when I reference it specifically?

 

 

T2[a, b, matrix]

T2[a, b] = (Matrix(4, 4, {(1, 1) = 1, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (2, 1) = 0, (2, 2) = -1, (2, 3) = 0, (2, 4) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = -1, (3, 4) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = -1}))

(16)

T2[3, 3]

-1

(17)

 

Lower indices work as expected.

 

``



Download TensorArray_Tetrads.mw

 

 

 

 

@ecterrab 

Thanks for the update. I'll keep updating and let you know if I find any issues or bugs. :)

@ecterrab 

Hi Edgardo, 

The physics package you are developing is amazing. I have looked through the null tetrad worksheet you suggested and it seems to contain everything I was looking for. 

I will stop hasseling you, and let you finish the package in peace, however I was wondering if there was a changelog or something I could subscribe to, in order to keep up to date with your progress. 

One final thing, I was unable to use the Tetrad sub-package, there seems to be a bug or perhaps something wrong with my setup. I'm using the linux version of Maple 18 (fully updated). And I have this:

> with(Physics): <- All good

>with(Physics:-Tetrads); - Error, (in with) user level initialization for package `Physics:-Tetrads' failed: Print is not a command in the Physics package

Oh, I should mention the physics package was updated as of today. 

Thanks again.

@ecterrab 

Hi Edgardo, 

I have seen the SumOverRepeateIndicies (I dont use this as I am summing over different ranges of objects). 

If I have a non-minkwoskian metric, the inverse metric g_[~mu,~nu] is obviously different to the metric g_[mu,nu]. 

I have tensors of this form, where adding the tilde makes life much easier as it automatically does the contraction with the metric for me. 

The command 

> add(g_[~mu,2], ~mu=1..2) does the addition g_[1,2] + g_[2,2] not the desired addition g_[~1,2] + g_[~2,2].  

There may not be an easy way of doing this, and I may be able to reformulate my problem all in terms of tensors. 

 

1 2 Page 1 of 2