Lucero has asked for the wisdom of the Perl Monks concerning the following question:
Hello. Im trying to calculate the transition probability for a constant pertubation.
What I have to do is apply the RK4 method to a vectorial function
:fn(cm, t) =sum over m((i/h)*<nVm>exp(i(Em En)t/h))Cm)
The values of cm^2 must oscillate around 1 and 0 but I run the code the values I get are bigger and dont have the oscillatory behavior.
use warnings;
use strict;
use PDL;
use PDL::Complex;
my $eprim=0.001;
my $X =pdl [
[ 0, sqrt(1),0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
+],
[sqrt(1),0, sqrt(2),0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
+,0,0],
[0, sqrt(2),0,sqrt(3),0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
+,0,0],
[0,0,sqrt(3),0,sqrt(4),0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+0,0],
[0,0,0,sqrt(4),0,sqrt(5),0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+0,0],
[0,0,0,0,sqrt(5),0,sqrt(6),0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+0,0],
[0,0,0,0,0,sqrt(6),0,sqrt(7),0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+0,0],
[0,0,0,0,0,0,sqrt(7),0,sqrt(8),0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+0,0],
[0,0,0,0,0,0,0,sqrt(8),0,sqrt(9),0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+0,0],
[0,0,0,0,0,0,0,0,sqrt(9),0,sqrt(10),0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
+,0,0],
[0,0,0,0,0,0,0,0,0,sqrt(10),0,sqrt(11),0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+0,0,0],
[0,0,0,0,0,0,0,0,0,0,sqrt(11),0,sqrt(12),0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,sqrt(12),0,sqrt(13),0,0,0,0,0,0,0,0,0,0,0,0,0,
+0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,sqrt(13),0,sqrt(14),0,0,0,0,0,0,0,0,0,0,0,0,
+0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,sqrt(14),0,sqrt(15),0,0,0,0,0,0,0,0,0,0,
+0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,sqrt(15),0,sqrt(16),0,0,0,0,0,0,0,0,0,
+0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,sqrt(16),0,sqrt(17),0,0,0,0,0,0,0,0,
+0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,sqrt(17),0,sqrt(18),0,0,0,0,0,0,0,
+0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,sqrt(18),0,sqrt(19),0,0,0,0,0,0,
+0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,sqrt(19),0,sqrt(20),0,0,0,0,0,
+0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,sqrt(20),0,sqrt(21),0,0,0,0,
+0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,sqrt(21),0,sqrt(22),0,0,0,
+0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,sqrt(22),0,sqrt(23),0,0,
+0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,sqrt(23),0,sqrt(24),0,
+0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,sqrt(24),0,sqrt(25),
+0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,sqrt(25),0,sqrt(26
+),0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,sqrt(26),0,sqrt(
+27),0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,sqrt(27),0,sqr
+t(27)],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,sqrt(28),0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,sqrt(28)]
];
my $vec= pdl [
[ 0.4037003675364034,0.0000000000000169,0.8381213708160592,0.0000
+000000000158,0.3190879932741590,0.0000000000000003,0.169355258154224
+8,0.0000000000000032,0.0251289794138005,0.0000000000000012,0.053630
+0581595280,0.0000000000000003,0.0217496956429717,0.0000000000000005,
+0.0021651419544516,0.0000000000000002,0.0080249333137898,0.00000000
+00000000,0.0051470455076793,0.0000000000000001,0.0013650508619321,0
+.0000000000000001,0.0006214095672646,0.0000000000000000,0.0010251417
+156844,0.0000000000000000,0.0007344546098787,0.0000000000000000,0.
+0003358368423501,0.0000000000000000 ],
[ 0.0000000000000095,0.7226441842331763,0.0000000000000196,0.674369
+4628487641,0.0000000000000074,0.0126709869133627,0.0000000000000039
+,0.1390097337149938,0.0000000000000004,0.0532637295249766,0.0000000
+000000012,0.0108058217041577,0.0000000000000004,0.0211856923535512,
+0.0000000000000004,0.0104565487807121,0.0000000000000001,0.000802132
+1669403,0.0000000000000003,0.0027244498764788,0.0000000000000004,0.
+0023881197288006,0.0000000000000006,0.0010161248442056,0.0000000000
+000009,0.0000173027596815,0.0000000000000008,0.0003416491521171,0.00
+00000000000006,0.0002612727358195 ],
[ 0.7265845261792991,0.0000000000000004,0.0822771629823105,0.000000
+0000000005,0.6492960190820365,0.0000000000000003,0.0423366396915736
+,0.0000000000000002,0.1758256491790640,0.0000000000000000,0.09794090
+67513942,0.0000000000000001,0.0028801939655116,0.0000000000000000,
+0.0292811602848338,0.0000000000000000,0.0218748117671870,0.000000000
+0000000,0.0068869463649583,0.0000000000000000,0.0019130055021425,0
+.0000000000000000,0.0039694711332578,0.0000000000000000,0.0028100428
+723103,0.0000000000000000,0.0012454407121507,0.0000000000000000,0.
+0003466238255135,0.0000000000000000 ],
[ 0.0000000000000006,0.6010733799547240,0.0000000000000001,0.593736
+8797404066,0.0000000000000004,0.4447140685185609,0.0000000000000001
+,0.2753043611569276,0.0000000000000001,0.0177297849332504,0.00000000
+00000001,0.0947825000175309,0.0000000000000001,0.0539543068519576,0
+.0000000000000003,0.0064734473245184,0.0000000000000001,0.013191845
+5031492,0.0000000000000004,0.0125678040256591,0.0000000000000000,0.0
+057077804477274,0.0000000000000002,0.0003357404158152,0.000000000000
+0007,0.0018880365733295,0.0000000000000005,0.0019443492967520,0.000
+0000000000001,0.0010789173843832 ],
[ 0.5430269760782945,0.0000000000000001,0.4474156032695855,0.00000
+00000000000,0.5801000179013242,0.0000000000000002,0.2404758383494775,
+0.0000000000000002,0.3088519136114667,0.0000000000000001,0.08237883
+44472166,0.0000000000000000,0.0541280464126355,0.0000000000000001,0
+.0660462396199596,0.0000000000000000,0.0305535579726857,0.0000000000
+000000,0.0009344967473301,0.0000000000000000,0.0099964129968316,0.
+0000000000000000,0.0087666379065681,0.0000000000000000,0.0040952318
+231305,0.0000000000000000,0.0007413411176992,0.0000000000000000,0.00
+03124230406161,0.0000000000000000 ],
[ 0.0000000000000000,0.3321617007160794,0.0000000000000001,0.331456
+4696269089,0.0000000000000002,0.7974495377138241,0.0000000000000004
+,0.0513734221888449,0.0000000000000001,0.3011617902587023,0.0000000
+000000004,0.2082370289532634,0.0000000000000003,0.0347605467840875,
+0.0000000000000001,0.0497057528447874,0.0000000000000006,0.05212963
+36393288,0.0000000000000008,0.0250539045206587,0.0000000000000008,0
+.0023068831550118,0.0000000000000003,0.0074112442261458,0.0000000000
+000001,0.0079651410760186,0.0000000000000001,0.0051288837146889,0.0
+000000000000001,0.0022324847785247 ],
[ 0.1172897082030818,0.0000000000000002,0.2964728523725000,0.000000
+0000000003,0.1841483714928239,0.0000000000000005,0.8047717594409407
+,0.0000000000000009,0.3242847252115117,0.0000000000000004,0.16364977
+62965390,0.0000000000000000,0.2529941067981705,0.0000000000000002,0.
+1294310921497701,0.0000000000000002,0.0062237278459524,0.0000000000
+000001,0.0444860334138650,0.0000000000000000,0.0406263992613470,0.0
+000000000000001,0.0185369919706947,0.0000000000000001,0.000753458698
+2924,0.0000000000000000,0.0063516252770288,0.0000000000000000,0.005
+0724227338856,0.0000000000000000 ],
[ 0.0000000000000000,0.0746816751641777,0.0000000000000003,0.279960
+3944858002,0.0000000000000004,0.0369884510514530,0.0000000000000016,
+0.7185111750625419,0.0000000000000002,0.5499552865074182,0.00000000
+00000006,0.0532991069401313,0.0000000000000002,0.2088524831915762,0.
+0000000000000011,0.1984245945881420,0.0000000000000011,0.0856956117
+134249,0.0000000000000004,0.0038508424287722,0.0000000000000007,0.03
+73833868060850,0.0000000000000012,0.0329127019021361,0.0000000000000
+010,0.0166509762582109,0.0000000000000004,0.0042911848981398,0.000
+0000000000001,0.0002821243295013 ],
[ 0.0214780579913734,0.0000000000000003,0.0372177058109391,0.000000
+0000000003,0.2944751450270211,0.0000000000000008,0.075683024140412
+2,0.0000000000000008,0.5412748833836341,0.0000000000000006,0.668676
+2866715919,0.0000000000000003,0.2932317876571968,0.0000000000000000,
+0.0735660926692254,0.0000000000000002,0.2050075536886243,0.00000000
+00000001,0.1583105998876948,0.0000000000000000,0.0602462844067277,
+0.0000000000000001,0.0115818910911991,0.0000000000000001,0.039716023
+0068139,0.0000000000000001,0.0369608788819537,0.0000000000000001,0.
+0202264690881708,0.0000000000000000 ],
[ 0.0000000000000000,0.0238139307663168,0.0000000000000000,0.0126527
+912037899,0.0000000000000002,0.3264304157504123,0.0000000000000002,
+0.1354442358079246,0.0000000000000006,0.3188188145305269,0.000000000
+0000004,0.6646875721923707,0.0000000000000002,0.4961929745611114,0.
+0000000000000003,0.1255681225264503,0.0000000000000002,0.127000931
+5120212,0.0000000000000003,0.1845765986762271,0.0000000000000001,0.
+1230738673852459,0.0000000000000004,0.0376175711869895,0.00000000000
+00006,0.0187532845837444,0.0000000000000005,0.0354892349606020,0.00
+00000000000004,0.0238711615056090 ],
[ 0.0017318840561273,0.0000000000000002,0.0324109249717848,0.000000
+0000000007,0.0699841577530595,0.0000000000000002,0.349302010000345
+1,0.0000000000000011,0.1121863328196896,0.0000000000000014,0.10921
+03091430736,0.0000000000000000,0.5420611831710406,0.0000000000000013
+,0.6097584088048894,0.0000000000000015,0.3508901747526678,0.0000000
+000000009,0.0471175681506413,0.0000000000000002,0.1258773221542367
+,0.0000000000000002,0.1564524321190749,0.0000000000000002,0.1120328
+335999530,0.0000000000000002,0.0570522771318334,0.0000000000000001,
+0.0203908274536104,0.0000000000000001 ],
[ 0.0000000000000000,0.0016213341083221,0.0000000000000001,0.051098
+7268286506,0.0000000000000003,0.1334540474130648,0.0000000000000009,
+0.3479817094912636,0.0000000000000007,0.0311921491179282,0.00000000
+00000001,0.0444673033625089,0.0000000000000014,0.3398228370041564,0.
+0000000000000012,0.5916539149462523,0.0000000000000005,0.5108305696
+513560,0.0000000000000002,0.2420842730477737,0.0000000000000003,0
+.0103981442743016,0.0000000000000013,0.1459764073717070,0.0000000000
+000019,0.1690313714662625,0.0000000000000019,0.1276772440329932,0.0
+000000000000014,0.0644813878617362 ],
[ 0.0011147979892237,0.0000000000000003,0.0124527283486975,0.0000000
+000000002,0.0947114844587487,0.0000000000000008,0.2372976408624008,0
+.0000000000000005,0.3557847566629774,0.0000000000000003,0.0635298085
+436199,0.0000000000000005,0.1713173200400193,0.0000000000000002,0.
+0972899291799430,0.0000000000000004,0.4717862901085890,0.00000000000
+00011,0.5720800153535275,0.0000000000000016,0.3811052943887194,0.00
+00000000000018,0.0912455512043619,0.0000000000000014,0.119297875587
+9619,0.0000000000000009,0.1817055145006404,0.0000000000000005,0.120
+2320521959933,0.0000000000000002 ],
[ 0.0000000000000000,0.0034343034398039,0.0000000000000001,0.0299
+555389414543,0.0000000000000004,0.1423655878081122,0.0000000000000
+012,0.2861228319386651,0.0000000000000016,0.2698629683232922,0.00
+00000000000006,0.1808943245053649,0.0000000000000014,0.13910518243619
+62,0.0000000000000000,0.0647414559720777,0.0000000000000010,0.2908
+420542400688,0.0000000000000007,0.5446044409698884,0.0000000000000011
+,0.5243579846618940,0.0000000000000024,0.3154050447298352,0.0000000
+000000033,0.0894563929271685,0.0000000000000028,0.0389334639461899
+,0.0000000000000014,0.0522795627061793 ],
[ 0.0002173011090244,0.0000000000000015,0.0060764847307131,0.0000
+000000000005,0.0347032140814115,0.0000000000000003,0.12141153092802
+55,0.0000000000000004,0.1627557440827222,0.0000000000000027,0.08257
+98997978443,0.0000000000000007,0.1886866183315523,0.000000000000000
+9,0.0554160985866513,0.0000000000000009,0.0016082286660558,0.00000
+00000000035,0.1445886337968438,0.0000000000000039,0.373441057958475
+9,0.0000000000000024,0.5230941284069114,0.0000000000000000,0.5241670
+456105690,0.0000000000000016,0.3997788067656795,0.0000000000000018,
+0.2084878025028327,0.0000000000000011 ],
[ 0.0000000000000000,0.0018740215138967,0.0000000000000000,0.020914
+3122794822,0.0000000000000002,0.0937250582686831,0.0000000000000006
+,0.2263227689510601,0.0000000000000010,0.2136068393701215,0.0000000
+000000003,0.0343030148937000,0.0000000000000006,0.2374515116591406,
+0.0000000000000009,0.1543298788071321,0.0000000000000014,0.048410648
+8830808,0.0000000000000011,0.0556695757276321,0.0000000000000009,0
+.1585122753165241,0.0000000000000015,0.4164529081832608,0.0000000000
+000016,0.5455319030422324,0.0000000000000021,0.4850082730133928,0.0
+000000000000015,0.2751494973249398 ],
[ 0.0003983422749079,0.0000000000000003,0.0084162036010590,0.000000
+0000000008,0.0612865660082999,0.0000000000000003,0.2247020891769621,0
+.0000000000000001,0.4305239272210992,0.0000000000000002,0.3202309758
+153931,0.0000000000000001,0.1918139117119362,0.0000000000000002,0.
+3095871293039516,0.0000000000000002,0.2464518281395873,0.00000000000
+00006,0.1656482308020451,0.0000000000000006,0.3359183916702893,0.00
+00000000000000,0.1094739661795629,0.0000000000000012,0.24302905217495
+86,0.0000000000000017,0.4099538699310810,0.0000000000000015,0.29157
+77578643249,0.0000000000000009 ],
[ 0.0000000000000000,0.0014954639775851,0.0000000000000000,0.01796859
+65171856,0.0000000000000000,0.0988989956441026,0.0000000000000001,0.2
+952876521886211,0.0000000000000003,0.4627556007193223,0.0000000000000
+004,0.2320927274343102,0.0000000000000000,0.2873578852921854,0.0000
+000000000001,0.2359945232765174,0.0000000000000001,0.311680124031071
+2,0.0000000000000003,0.0798333562246646,0.0000000000000004,0.341728
+8746164391,0.0000000000000003,0.1906004227513239,0.0000000000000001,
+0.1610665397841304,0.0000000000000002,0.3676708309232274,0.00000000
+00000002,0.2804853923818191 ],
[ 0.0001084747877321,0.0000000000000000,0.0027208960217152,0.00000000
+00000002,0.0242704434475306,0.0000000000000001,0.1153670170379691,0
+.0000000000000005,0.3193789054222693,0.0000000000000001,0.4908523839
+275092,0.0000000000000001,0.2761691414257316,0.0000000000000000,0.26
+00143419201901,0.0000000000000001,0.3068704363306086,0.000000000000
+0002,0.2542843527027897,0.0000000000000000,0.2086903757180524,0.0000
+000000000002,0.3365315876512961,0.0000000000000002,0.01653724898975
+50,0.0000000000000002,0.3099014215132993,0.0000000000000001,0.29551
+75636216057,0.0000000000000000 ],
[ 0.0000000000000000,0.0003844645253044,0.0000000000000000,0.0055
+700399702677,0.0000000000000000,0.0384100619532071,0.0000000000000
+000,0.1533025985895169,0.0000000000000000,0.3681698917672083,0.00
+00000000000001,0.4889286723129618,0.0000000000000002,0.19429626416
+19352,0.0000000000000001,0.3158440934858354,0.0000000000000001,0.24
+96900204706467,0.0000000000000001,0.2959512116085838,0.00000000000
+00000,0.1552508901502876,0.0000000000000001,0.3427579648786467,0.0
+000000000000002,0.0541435546335990,0.0000000000000000,0.28760313265
+88631,0.0000000000000001,0.2927760031123848 ],
[ 0.0000181563100889,0.0000000000000002,0.0005488425504935,0.0000000
+000000000,0.0060688055327925,0.0000000000000002,0.0372631786179999,
+0.0000000000000002,0.1425926143202678,0.0000000000000000,0.347871001
+1662555,0.0000000000000002,0.5076314191518138,0.0000000000000006,0.31
+23574515561076,0.0000000000000001,0.2091284023204721,0.000000000000
+0002,0.3689178831666572,0.0000000000000000,0.1344053015043881,0.000
+0000000000001,0.3299092963194629,0.0000000000000004,0.20822789510253
+15,0.0000000000000003,0.2240489250735314,0.0000000000000001,0.3188
+632967650468,0.0000000000000002 ],
[ 0.0000000000000000,0.0000646482239263,0.0000000000000000,0.001138
+5375166878,0.0000000000000000,0.0098338922987181,0.0000000000000000
+,0.0513935599025939,0.0000000000000001,0.1741934371091206,0.0000000
+000000002,0.3836429445985159,0.0000000000000004,0.5045604175490942,
+0.0000000000000004,0.2522551752686503,0.0000000000000001,0.257944102
+5360375,0.0000000000000004,0.3397702320139469,0.0000000000000003,0.1
+752821352580565,0.0000000000000001,0.3071408379256653,0.000000000000
+0001,0.2279977492913131,0.0000000000000002,0.2051227648591225,0.00
+00000000000004,0.3168989287342151 ],
[ 0.0000019661699240,0.0000000000000002,0.0000719533317298,0.00000000
+00000001,0.0009846334801556,0.0000000000000002,0.0077074332169402,0
+.0000000000000004,0.0392493847315142,0.0000000000000002,0.13658602674
+96776,0.0000000000000000,0.3260488848298075,0.0000000000000001,0.507
+5095625005354,0.0000000000000000,0.4189127624432365,0.00000000000000
+03,0.0412751650984030,0.0000000000000002,0.4061413387987063,0.0000
+000000000002,0.1232820067450982,0.0000000000000001,0.352583981765542
+7,0.0000000000000002,0.1125603524749182,0.0000000000000002,0.34863
+08267332468,0.0000000000000004 ],
[ 0.0000000000000000,0.0000071522511841,0.0000000000000000,0.0001
+531738843782,0.0000000000000000,0.0016457996800569,0.0000000000000
+000,0.0110351632854802,0.0000000000000000,0.0502042359369793,0.00
+00000000000001,0.1598827125229744,0.0000000000000002,0.35410559958
+16644,0.0000000000000003,0.5131002961113091,0.0000000000000004,0.
+3827211429896542,0.0000000000000001,0.0868953534697583,0.00000000000
+00002,0.4060663544932935,0.0000000000000002,0.0906849061585278,0.0000
+000000000001,0.3565218256839784,0.0000000000000002,0.09547988946696
+47,0.0000000000000000,0.3460718461864576 ],
[ 0.0000001292972631,0.0000000000000001,0.0000057407224638,0.00000
+00000000000,0.0000969677496746,0.0000000000000000,0.000957639115608
+1,0.0000000000000002,0.0063336604552633,0.0000000000000000,0.02981
+94086989161,0.0000000000000001,0.1025898937974792,0.0000000000000004
+,0.2581306279126550,0.0000000000000004,0.4608777058952639,0.0000000
+000000003,0.5289025115880573,0.0000000000000000,0.2501979027713430,
+0.0000000000000003,0.2459229288491326,0.0000000000000000,0.39768907
+21949536,0.0000000000000001,0.0441899989974140,0.0000000000000001,
+0.3828591066315501,0.0000000000000001 ],
[ 0.0000000000000000,0.0000004836289753,0.0000000000000000,0.00001259
+35691119,0.0000000000000000,0.0001674328686760,0.0000000000000000,0.0
+014204630187498,0.0000000000000000,0.0084233556059741,0.0000000000000
+000,0.0364726872642141,0.0000000000000000,0.1173353239240527,0.000000
+0000000001,0.2791623093846674,0.0000000000000001,0.4743110004921181,0
+.0000000000000002,0.5168194754141317,0.0000000000000002,0.21722128100
+95154,0.0000000000000001,0.2660675030932704,0.0000000000000000,0.38
+68615369372250,0.0000000000000001,0.0582093189941330,0.0000000000000
+002,0.3792058820414961 ],
[ 0.0000000043467313,0.0000000000000002,0.0000002354338414,0.0000000
+000000001,0.0000049191205605,0.0000000000000001,0.0000611103982443,0
+.0000000000000001,0.0005191484660093,0.0000000000000002,0.00322513771
+93584,0.0000000000000002,0.0151840473609755,0.0000000000000001,0.05
+51258405915315,0.0000000000000001,0.1547931322781031,0.00000000000000
+01,0.3318050095086971,0.0000000000000004,0.5216362528850688,0.00000
+00000000001,0.5386522451902931,0.0000000000000004,0.2190409351878887
+,0.0000000000000002,0.2724800121518238,0.0000000000000002,0.4221250
+094454858,0.0000000000000001 ],
[ 0.0000000000000000,0.0000000167699136,0.0000000000000000,0.0000
+005331602679,0.0000000000000000,0.0000087760999773,0.0000000000000
+000,0.0000937543740401,0.0000000000000000,0.0007150017363018,0.00
+00000000000000,0.0040915831845336,0.0000000000000000,0.01805307191
+06190,0.0000000000000000,0.0621905012332527,0.0000000000000000,0.
+1672502296997740,0.0000000000000001,0.3457819876105404,0.000000000
+0000002,0.5268406404622072,0.0000000000000002,0.5273801655434678,
+0.0000000000000002,0.1998433947965230,0.0000000000000000,0.280500487
+4860440,0.0000000000000003,0.4170641909080779 ],
[ 0.0000000000469012,0.0000000000000001,0.0000000031569888,0.0000000
+000000001,0.0000000829508766,0.0000000000000000,0.0000013139183774,0
+.0000000000000001,0.0000144662736142,0.0000000000000001,0.00011879541
+44757,0.0000000000000001,0.0007576595856397,0.0000000000000002,0.003
+8457059058298,0.0000000000000001,0.0157544770583618,0.00000000000000
+02,0.0524034914224598,0.0000000000000002,0.1412730634268711,0.000000
+0000000000,0.3051191205191895,0.0000000000000000,0.5135323458897670,
+0.0000000000000002,0.6319248579965772,0.0000000000000001,0.469989404
+2726668,0.0000000000000000 ],
[ 0.0000000000000000,0.0000000001872384,0.0000000000000000,0.000000
+0073974097,0.0000000000000000,0.0000001531161950,0.0000000000000000
+,0.0000020854332860,0.0000000000000000,0.0000206106783438,0.0000000
+000000000,0.0001558996179038,0.0000000000000000,0.0009318666076831,
+0.0000000000000000,0.0044889011562547,0.0000000000000000,0.017620794
+6824076,0.0000000000000000,0.0565953635430333,0.0000000000000000,0.
+1482728021860708,0.0000000000000000,0.3129386149972143,0.0000000000
+000000,0.5172629611893457,0.0000000000000001,0.6280792712999990,0.0
+000000000000001,0.4631473700646182 ]
]; #matrix T de los vectores del doblepozo
my $vect= transpose($vec); #matrix con los vectores en la columna
# <nVm>
my $V= $vec x $X x $vect;
my $step=0.1;
my $hbar=1;
#print &Vnm(0,1);
#print&Vnm(1,0);
my @Em=(3.49090786,3.48609059,0.94249804,0.72223315,0.84486159,1.9
+8334101,3.52168993,5.19562490,7.00900149,8.96758961,11.06144444445427
+60,13.1282317263938650,15.4597663541535190,18.2923955058715430,19.575
+1133106800170,21.3254642902799140,27.3060148356095880,33.615327480770
+2470,48.2864949409823790,58.1657460989797170,83.0782032408711330,97.3
+475523518670660,137.5307666823436800,157.3985331875067300,221.5161226
+004345300,248.6116702653141900,353.3740736965299900,390.0680354316574
+500,577.9792266930371600,628.4681868621187300);
my $w=$Em[1]$Em[0];
sub Vnm {
my ($x1,$x2)=@_;
return $V>range([$x1,$x2]);
};
my $n;
my $m;
my $f;
my @c = (1, 0, 0, 0); #INITIAL CONDITIONS
#@c=F=(F0,F1,F2,F3)
sub F0 {
my ($t,@c)=@_;
$m=0;
$f=0;
for ($n=0;$n<4;$n++){
my $a= &Vnm($n,$m)*exp((i)*($Em[$m]$Em[$n])*($t/$hbar));
if ($n!=$m){
$f = $f + $a*((i)*($eprim)*cos($w*$t)*$c[$n]);
};
};
return $f;
};
sub F1 {
my ($t,@c)=@_;
$m=1;
$f=0;
for ($n=0;$n<4;$n++){
my $a= &Vnm($n,$m)*exp((i)*($Em[$m]$Em[$n])*($t/$hbar));
if ($n!=$m){
$f = $f + $a*((i)*($eprim)*cos($w*$t)*$c[$n]);
};
};
return $f;
};
sub F2 {
my ($t,@c)=@_;
$m=2;
$f=0;
for ($n=0;$n<4;$n++){
my $a= &Vnm($n,$m)*exp((i)*($Em[$m]$Em[$n])*($t/$hbar));
if ($n!=$m){
$f = $f + $a*(i)*($eprim)*cos($w*$t)*$c[$n];
};
};
return $f;
};
sub F3 {
my ($t,@c)=@_;
$m=3;
$f=0;
for ($n=0;$n<4;$n++){
my $a= &Vnm($n,$m)*exp((i)*($Em[$m]$Em[$n])*($t/$hbar));
for ($n!=$m){
$f = $f +$a*((i)*($eprim)*cos($w*$t)*$c[$n]);
};
};
return $f;
};
#Método RK4
my $k1;
my $k2;
my $k3;
my $k4;
my $k1_1;
my $k2_1;
my $k3_1;
my $k4_1;
my $k1_2;
my $k2_2;
my $k3_2;
my $k4_2;
my $k1_3;
my $k2_3;
my $k3_3;
my $k4_3;
sub RK4_c0 {
my ($t,@c)=@_;
sub k1 { my ($t1,@c)=@_;return $step*&F0($t1,@c); };
sub k2 { my ($t2,@c)=@_;return $step*&F0(($t2+$step/2),( @c + (&k1($t2
+)/2)));};
sub k3 { my ($t3,@c)=@_;return $step*&F0(($t3+$step/2),( @c + (&k2($t3
+)/2))) ;};
sub k4 { my ($t4,@c)=@_;return $step*&F0(($t4+$step),( @c + &k3($t4)))
+;};
return $c[0] + (1/6)*(&k1($t,@c)+ 2*(&k2($t,@c) + &k3($t,@c)) + &k4($
+t,@c));
};
sub RK4_c1 {
my ($t,@c)=@_;
sub k1_1 { my ($t1,@c)=@_;return $step*&F1($t1,@c); };
sub k2_1 { my ($t2,@c)=@_;return $step*&F1(($t2+$step/2),( @c + &k1($t
+2/2)));};
sub k3_1 { my ($t3,@c)=@_;return $step*&F1(($t3+$step/2),( @c + &k2($t
+3)/2)) ;};
sub k4_1 { my ($t4,@c)=@_;return $step*&F1(($t4+$step),( @c + &k3($t4)
+));};
return $c[1] + (1/6)*(&k1_1($t,@c)+ 2*(&k2_1($t,@c) + &k3_1($t,@c)) +
+ &k4_1($t,@c));
};
sub RK4_c2 {
my ($t,@c)=@_;
sub k1_2 { my ($t1,@c)=@_;return $step*&F2($t1,@c); };
sub k2_2 { my ($t2,@c)=@_;return $step*&F2(($t2+$step/2),( @c + &k1($t
+2)/2));};
sub k3_2 { my ($t3,@c)=@_;return $step*&F2(($t3+$step/2),( @c + &k2($t
+3)/2)) ;};
sub k4_2 { my ($t4,@c)=@_;return $step*&F2(($t4+$step),( @c + &k3($t4)
+));};
return $c[2] + (1/6)*(&k1_2($t,@c)+ 2*(&k2_2($t,@c) + &k3_2($t,@c)) +
+ &k4_2($t,@c));
};
sub RK4_c3 {
my ($t,@c)=@_;
sub k1_3 { my ($t1,@c)=@_;return $step*&F3($t1,@c); };
sub k2_3 { my ($t2,@c)=@_;return $step*&F3(($t2+$step/2),( @c + &k1($t
+2)/2));};
sub k3_3 { my ($t3,@c)=@_;return $step*&F3(($t3+$step/2),( @c + &k2($t
+3)/2)) ;};
sub k4_3 { my ($t4,@c)=@_;return $step*&F3(($t4+$step),( @c + &k3($t4)
+));};
return $c[3] + (1/6)*(&k1_3($t,@c)+ 2*(&k2_3($t,@c) + &k3_3($t,@c)) +
+ &k4_3($t,@c));
};
my $t=0;
open (FILE , ">C0_V1W01p2.dat");
for ($t=0; $t<2500;$t+=$step){
#print " @c $t \n";
#print abs($c[0])*abs($c[0]), "\n";
$c[0]= &RK4_c0($t,@c);
$c[1]= &RK4_c1($t,@c);
$c[2]= &RK4_c2($t,@c);
$c[3]= &RK4_c3($t,@c);
my $x= abs($c[0])**2;
my $b= abs($c[1])**2;
my $s= abs($c[2])**2;
my $d= abs($c[3])**2;
print FILE "$x $b $s $d $t\n";
};
close (FILE);
Re: RK4 for a vectorial function
by Fletch (Chancellor) on Nov 26, 2021 at 08:51 UTC

Another helping you get help suggestion: you might provide more context what an "RK4 for a vectorial function" is, because I for sure have no clue what you're referring to off hand. If this is some common algorithm for some specific problem domain that's all well and good, but the peanut gallery here is (like myself) probably completely unfamiliar with what you're talking about.
I mean a wild guess and wolfram alpha query turns up RungeKutta Method which looks similar if I squint and scroll around. But giving that context out up front rather than making people trying to help do the work first helps us (well, someone here maybe) help you.
The cake is a lie.
The cake is a lie.
The cake is a lie.
 [reply] 
Re: RK4 for a vectorial function
by vr (Curate) on Nov 26, 2021 at 07:12 UTC

I definitely understand nothing about what this all is supposed to do. At 1st glance: are you sure about calling
&F0($t1,@c);
...
&F0(($t2+$step/2),( @c + (&k1($t2)/2)));
? What (how many) arguments do you think F0 receives in each case?
Further, k2_1 calls k1. Just for sake of symmetry, I wonder if it should have been k1_1. Have you tested/debugged each subroutine/section, so they work as expected, instead of typing this huge and complex code to observe final results? And I suspect other monks' advice about code and PBP will follow.
 [reply] [d/l] [select] 
Re: RK4 for a vectorial function
by jwkrahn (Monsignor) on Nov 27, 2021 at 02:24 UTC

#!/usr/bin/perl
use warnings;
use strict;
use PDL;
use PDL::Complex;
my $eprim = 0.001;
my $X = pdl( [
[ 0, sqrt( 1 ), 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0,
+0 ],
[ sqrt( 1 ), 0, sqrt( 2 ), 0, 0, 0,
+ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0,
+0 ],
[ 0, sqrt( 2 ), 0, sqrt( 3 ), 0, 0,
+ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0,
+0 ],
[ 0, 0, sqrt( 3 ), 0, sqrt( 4 ), 0,
+ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0,
+0 ],
[ 0, 0, 0, sqrt( 4 ), 0, sqrt( 5
+), 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0,
+0 ],
[ 0, 0, 0, 0, sqrt( 5 ), 0,
+ sqrt( 6 ), 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0,
+0 ],
[ 0, 0, 0, 0, 0, sqrt( 6
+), 0, sqrt( 7 ), 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0,
+0 ],
[ 0, 0, 0, 0, 0, 0,
+ sqrt( 7 ), 0, sqrt( 8 ), 0, 0, 0,
+ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0,
+0 ],
[ 0, 0, 0, 0, 0, 0,
+ 0, sqrt( 8 ), 0, sqrt( 9 ), 0, 0,
+ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0,
+0 ],
[ 0, 0, 0, 0, 0, 0,
+ 0, 0, sqrt( 9 ), 0, sqrt( 10 ), 0,
+ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0,
+0 ],
[ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, sqrt( 10 ), 0, sqrt( 11
+), 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0,
+0 ],
[ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, sqrt( 11 ), 0,
+ sqrt( 12 ), 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0,
+0 ],
[ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, sqrt( 12
+), 0, sqrt( 13 ), 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0,
+0 ],
[ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0,
+ sqrt( 13 ), 0, sqrt( 14 ), 0, 0, 0,
+ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0,
+0 ],
[ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0,
+ 0, 0, sqrt( 14 ), 0, sqrt( 15 ), 0,
+ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0,
+0 ],
[ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, sqrt( 15 ), 0, sqrt(
+16 ), 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0,
+0 ],
[ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, sqrt( 16 ), 0,
+ sqrt( 17 ), 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0,
+0 ],
[ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, sqrt(
+17 ), 0, sqrt( 18 ), 0, 0, 0, 0,
+ 0, 0, 0, 0, 0,
+0 ],
[ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0,
+ sqrt( 18 ), 0, sqrt( 19 ), 0, 0, 0,
+ 0, 0, 0, 0, 0,
+0 ],
[ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0,
+ 0, sqrt( 19 ), 0, sqrt( 20 ), 0, 0,
+ 0, 0, 0, 0, 0,
+0 ],
[ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0,
+ 0, 0, sqrt( 20 ), 0, sqrt( 21 ), 0,
+ 0, 0, 0, 0, 0,
+0 ],
[ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, sqrt( 21 ), 0, sqr
+t( 22 ), 0, 0, 0, 0, 0,
+0 ],
[ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, sqrt( 22 ), 0,
+ sqrt( 23 ), 0, 0, 0, 0,
+0 ],
[ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, sqr
+t( 23 ), 0, sqrt( 24 ), 0, 0, 0,
+0 ],
[ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0,
+ sqrt( 24 ), 0, sqrt( 25 ), 0, 0,
+0 ],
[ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0,
+ 0, sqrt( 25 ), 0, sqrt( 26 ), 0,
+0 ],
[ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0,
+ 0, 0, sqrt( 26 ), 0, sqrt( 27 ),
+0 ],
[ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, sqrt( 27 ), 0,
+sqrt( 27 ) ],
[ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, sqrt( 28 ),
+0 ],
[ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0,
+sqrt( 28 ) ],
] );
my $vec = pdl( [
[ 0.4037003675364034, 0.0000000000000169, 0.8381213708160592, 
+0.0000000000000158, 0.3190879932741590, 0.0000000000000003, 0.1693
+552581542248, 0.0000000000000032, 0.0251289794138005, 0.0000000000
+000012, 0.0536300581595280, 0.0000000000000003, 0.0217496956429717
+, 0.0000000000000005, 0.0021651419544516, 0.0000000000000002, 0.0
+080249333137898, 0.0000000000000000, 0.0051470455076793, 0.0000000
+000000001, 0.0013650508619321, 0.0000000000000001, 0.0006214095672
+646, 0.0000000000000000, 0.0010251417156844, 0.0000000000000000, 
+0.0007344546098787, 0.0000000000000000, 0.0003358368423501, 0.0000
+000000000000 ],
[ 0.0000000000000095, 0.7226441842331763, 0.0000000000000196,
+0.6743694628487641, 0.0000000000000074, 0.0126709869133627, 0.0000
+000000000039, 0.1390097337149938, 0.0000000000000004, 0.0532637295
+249766, 0.0000000000000012, 0.0108058217041577, 0.0000000000000004
+, 0.0211856923535512, 0.0000000000000004, 0.0104565487807121, 0.0
+000000000000001, 0.0008021321669403, 0.0000000000000003, 0.0027244
+498764788, 0.0000000000000004, 0.0023881197288006, 0.0000000000000
+006, 0.0010161248442056, 0.0000000000000009, 0.0000173027596815, 
+0.0000000000000008, 0.0003416491521171, 0.0000000000000006, 0.0002
+612727358195 ],
[ 0.7265845261792991, 0.0000000000000004, 0.0822771629823105, 
+0.0000000000000005, 0.6492960190820365, 0.0000000000000003, 0.0423
+366396915736, 0.0000000000000002, 0.1758256491790640, 0.0000000000
+000000, 0.0979409067513942, 0.0000000000000001, 0.0028801939655116
+, 0.0000000000000000, 0.0292811602848338, 0.0000000000000000, 0.0
+218748117671870, 0.0000000000000000, 0.0068869463649583, 0.0000000
+000000000, 0.0019130055021425, 0.0000000000000000, 0.0039694711332
+578, 0.0000000000000000, 0.0028100428723103, 0.0000000000000000,
+0.0012454407121507, 0.0000000000000000, 0.0003466238255135, 0.0000
+000000000000 ],
[ 0.0000000000000006, 0.6010733799547240, 0.0000000000000001, 
+0.5937368797404066, 0.0000000000000004, 0.4447140685185609, 0.0000
+000000000001, 0.2753043611569276, 0.0000000000000001, 0.0177297849
+332504, 0.0000000000000001, 0.0947825000175309, 0.0000000000000001
+, 0.0539543068519576, 0.0000000000000003, 0.0064734473245184, 0.0
+000000000000001, 0.0131918455031492, 0.0000000000000004, 0.0125678
+040256591, 0.0000000000000000, 0.0057077804477274, 0.0000000000000
+002, 0.0003357404158152, 0.0000000000000007, 0.0018880365733295, 
+0.0000000000000005, 0.0019443492967520, 0.0000000000000001, 0.0010
+789173843832 ],
[ 0.5430269760782945, 0.0000000000000001, 0.4474156032695855, 
+0.0000000000000000, 0.5801000179013242, 0.0000000000000002, 0.2404
+758383494775, 0.0000000000000002, 0.3088519136114667, 0.0000000000
+000001, 0.0823788344472166, 0.0000000000000000, 0.0541280464126355
+, 0.0000000000000001, 0.0660462396199596, 0.0000000000000000, 0.0
+305535579726857, 0.0000000000000000, 0.0009344967473301, 0.0000000
+000000000, 0.0099964129968316, 0.0000000000000000, 0.0087666379065
+681, 0.0000000000000000, 0.0040952318231305, 0.0000000000000000,
+0.0007413411176992, 0.0000000000000000, 0.0003124230406161, 0.0000
+000000000000 ],
[ 0.0000000000000000, 0.3321617007160794, 0.0000000000000001, 
+0.3314564696269089, 0.0000000000000002, 0.7974495377138241, 0.0000
+000000000004, 0.0513734221888449, 0.0000000000000001, 0.3011617902
+587023, 0.0000000000000004, 0.2082370289532634, 0.0000000000000003
+, 0.0347605467840875, 0.0000000000000001, 0.0497057528447874, 0.0
+000000000000006, 0.0521296336393288, 0.0000000000000008, 0.0250539
+045206587, 0.0000000000000008, 0.0023068831550118, 0.0000000000000
+003, 0.0074112442261458, 0.0000000000000001, 0.0079651410760186,
+0.0000000000000001, 0.0051288837146889, 0.0000000000000001, 0.0022
+324847785247 ],
[ 0.1172897082030818, 0.0000000000000002, 0.2964728523725000,
+0.0000000000000003, 0.1841483714928239, 0.0000000000000005, 0.8047
+717594409407, 0.0000000000000009, 0.3242847252115117, 0.0000000000
+000004, 0.1636497762965390, 0.0000000000000000, 0.2529941067981705
+, 0.0000000000000002, 0.1294310921497701, 0.0000000000000002, 0.0
+062237278459524, 0.0000000000000001, 0.0444860334138650, 0.0000000
+000000000, 0.0406263992613470, 0.0000000000000001, 0.0185369919706
+947, 0.0000000000000001, 0.0007534586982924, 0.0000000000000000,
+0.0063516252770288, 0.0000000000000000, 0.0050724227338856, 0.0000
+000000000000 ],
[ 0.0000000000000000, 0.0746816751641777, 0.0000000000000003,
+0.2799603944858002, 0.0000000000000004, 0.0369884510514530, 0.0000
+000000000016, 0.7185111750625419, 0.0000000000000002, 0.5499552865
+074182, 0.0000000000000006, 0.0532991069401313, 0.0000000000000002
+, 0.2088524831915762, 0.0000000000000011, 0.1984245945881420, 0.0
+000000000000011, 0.0856956117134249, 0.0000000000000004, 0.0038508
+424287722, 0.0000000000000007, 0.0373833868060850, 0.0000000000000
+012, 0.0329127019021361, 0.0000000000000010, 0.0166509762582109, 
+0.0000000000000004, 0.0042911848981398, 0.0000000000000001, 0.0002
+821243295013 ],
[ 0.0214780579913734, 0.0000000000000003, 0.0372177058109391,
+0.0000000000000003, 0.2944751450270211, 0.0000000000000008, 0.0756
+830241404122, 0.0000000000000008, 0.5412748833836341, 0.0000000000
+000006, 0.6686762866715919, 0.0000000000000003, 0.2932317876571968
+, 0.0000000000000000, 0.0735660926692254, 0.0000000000000002, 0.2
+050075536886243, 0.0000000000000001, 0.1583105998876948, 0.0000000
+000000000, 0.0602462844067277, 0.0000000000000001, 0.0115818910911
+991, 0.0000000000000001, 0.0397160230068139, 0.0000000000000001,
+0.0369608788819537, 0.0000000000000001, 0.0202264690881708, 0.0000
+000000000000 ],
[ 0.0000000000000000, 0.0238139307663168, 0.0000000000000000,
+0.0126527912037899, 0.0000000000000002, 0.3264304157504123, 0.0000
+000000000002, 0.1354442358079246, 0.0000000000000006, 0.3188188145
+305269, 0.0000000000000004, 0.6646875721923707, 0.0000000000000002
+, 0.4961929745611114, 0.0000000000000003, 0.1255681225264503, 0.0
+000000000000002, 0.1270009315120212, 0.0000000000000003, 0.1845765
+986762271, 0.0000000000000001, 0.1230738673852459, 0.0000000000000
+004, 0.0376175711869895, 0.0000000000000006, 0.0187532845837444, 
+0.0000000000000005, 0.0354892349606020, 0.0000000000000004, 0.0238
+711615056090 ],
[ 0.0017318840561273, 0.0000000000000002, 0.0324109249717848,
+0.0000000000000007, 0.0699841577530595, 0.0000000000000002, 0.3493
+020100003451, 0.0000000000000011, 0.1121863328196896, 0.0000000000
+000014, 0.1092103091430736, 0.0000000000000000, 0.5420611831710406
+, 0.0000000000000013, 0.6097584088048894, 0.0000000000000015, 0.3
+508901747526678, 0.0000000000000009, 0.0471175681506413, 0.0000000
+000000002, 0.1258773221542367, 0.0000000000000002, 0.1564524321190
+749, 0.0000000000000002, 0.1120328335999530, 0.0000000000000002,
+0.0570522771318334, 0.0000000000000001, 0.0203908274536104, 0.0000
+000000000001 ],
[ 0.0000000000000000, 0.0016213341083221, 0.0000000000000001, 
+0.0510987268286506, 0.0000000000000003, 0.1334540474130648, 0.0000
+000000000009, 0.3479817094912636, 0.0000000000000007, 0.0311921491
+179282, 0.0000000000000001, 0.0444673033625089, 0.0000000000000014
+, 0.3398228370041564, 0.0000000000000012, 0.5916539149462523, 0.0
+000000000000005, 0.5108305696513560, 0.0000000000000002, 0.2420842
+730477737, 0.0000000000000003, 0.0103981442743016, 0.0000000000000
+013, 0.1459764073717070, 0.0000000000000019, 0.1690313714662625,
+0.0000000000000019, 0.1276772440329932, 0.0000000000000014, 0.0644
+813878617362 ],
[ 0.0011147979892237, 0.0000000000000003, 0.0124527283486975,
+0.0000000000000002, 0.0947114844587487, 0.0000000000000008, 0.2372
+976408624008, 0.0000000000000005, 0.3557847566629774, 0.0000000000
+000003, 0.0635298085436199, 0.0000000000000005, 0.1713173200400193
+, 0.0000000000000002, 0.0972899291799430, 0.0000000000000004, 0.4
+717862901085890, 0.0000000000000011, 0.5720800153535275, 0.0000000
+000000016, 0.3811052943887194, 0.0000000000000018, 0.0912455512043
+619, 0.0000000000000014, 0.1192978755879619, 0.0000000000000009,
+0.1817055145006404, 0.0000000000000005, 0.1202320521959933, 0.0000
+000000000002 ],
[ 0.0000000000000000, 0.0034343034398039, 0.0000000000000001, 
+0.0299555389414543, 0.0000000000000004, 0.1423655878081122, 0.0000
+000000000012, 0.2861228319386651, 0.0000000000000016, 0.2698629683
+232922, 0.0000000000000006, 0.1808943245053649, 0.0000000000000014
+, 0.1391051824361962, 0.0000000000000000, 0.0647414559720777, 0.0
+000000000000010, 0.2908420542400688, 0.0000000000000007, 0.5446044
+409698884, 0.0000000000000011, 0.5243579846618940, 0.0000000000000
+024, 0.3154050447298352, 0.0000000000000033, 0.0894563929271685, 
+0.0000000000000028, 0.0389334639461899, 0.0000000000000014, 0.0522
+795627061793 ],
[ 0.0002173011090244, 0.0000000000000015, 0.0060764847307131, 
+0.0000000000000005, 0.0347032140814115, 0.0000000000000003, 0.1214
+115309280255, 0.0000000000000004, 0.1627557440827222, 0.0000000000
+000027, 0.0825798997978443, 0.0000000000000007, 0.1886866183315523
+, 0.0000000000000009, 0.0554160985866513, 0.0000000000000009, 0.0
+016082286660558, 0.0000000000000035, 0.1445886337968438, 0.0000000
+000000039, 0.3734410579584759, 0.0000000000000024, 0.5230941284069
+114, 0.0000000000000000, 0.5241670456105690, 0.0000000000000016, 
+0.3997788067656795, 0.0000000000000018, 0.2084878025028327, 0.0000
+000000000011 ],
[ 0.0000000000000000, 0.0018740215138967, 0.0000000000000000,
+0.0209143122794822, 0.0000000000000002, 0.0937250582686831, 0.0000
+000000000006, 0.2263227689510601, 0.0000000000000010, 0.2136068393
+701215, 0.0000000000000003, 0.0343030148937000, 0.0000000000000006
+, 0.2374515116591406, 0.0000000000000009, 0.1543298788071321, 0.0
+000000000000014, 0.0484106488830808, 0.0000000000000011, 0.0556695
+757276321, 0.0000000000000009, 0.1585122753165241, 0.0000000000000
+015, 0.4164529081832608, 0.0000000000000016, 0.5455319030422324,
+0.0000000000000021, 0.4850082730133928, 0.0000000000000015, 0.2751
+494973249398 ],
[ 0.0003983422749079, 0.0000000000000003, 0.0084162036010590, 
+0.0000000000000008, 0.0612865660082999, 0.0000000000000003, 0.2247
+020891769621, 0.0000000000000001, 0.4305239272210992, 0.0000000000
+000002, 0.3202309758153931, 0.0000000000000001, 0.1918139117119362
+, 0.0000000000000002, 0.3095871293039516, 0.0000000000000002, 0.2
+464518281395873, 0.0000000000000006, 0.1656482308020451, 0.0000000
+000000006, 0.3359183916702893, 0.0000000000000000, 0.1094739661795
+629, 0.0000000000000012, 0.2430290521749586, 0.0000000000000017, 
+0.4099538699310810, 0.0000000000000015, 0.2915777578643249, 0.0000
+000000000009 ],
[ 0.0000000000000000, 0.0014954639775851, 0.0000000000000000,
+0.0179685965171856, 0.0000000000000000, 0.0988989956441026, 0.0000
+000000000001, 0.2952876521886211, 0.0000000000000003, 0.4627556007
+193223, 0.0000000000000004, 0.2320927274343102, 0.0000000000000000
+, 0.2873578852921854, 0.0000000000000001, 0.2359945232765174, 0.0
+000000000000001, 0.3116801240310712, 0.0000000000000003, 0.0798333
+562246646, 0.0000000000000004, 0.3417288746164391, 0.0000000000000
+003, 0.1906004227513239, 0.0000000000000001, 0.1610665397841304, 
+0.0000000000000002, 0.3676708309232274, 0.0000000000000002, 0.2804
+853923818191 ],
[ 0.0001084747877321, 0.0000000000000000, 0.0027208960217152,
+0.0000000000000002, 0.0242704434475306, 0.0000000000000001, 0.1153
+670170379691, 0.0000000000000005, 0.3193789054222693, 0.0000000000
+000001, 0.4908523839275092, 0.0000000000000001, 0.2761691414257316
+, 0.0000000000000000, 0.2600143419201901, 0.0000000000000001, 0.3
+068704363306086, 0.0000000000000002, 0.2542843527027897, 0.0000000
+000000000, 0.2086903757180524, 0.0000000000000002, 0.3365315876512
+961, 0.0000000000000002, 0.0165372489897550, 0.0000000000000002,
+0.3099014215132993, 0.0000000000000001, 0.2955175636216057, 0.0000
+000000000000 ],
[ 0.0000000000000000, 0.0003844645253044, 0.0000000000000000, 
+0.0055700399702677, 0.0000000000000000, 0.0384100619532071, 0.0000
+000000000000, 0.1533025985895169, 0.0000000000000000, 0.3681698917
+672083, 0.0000000000000001, 0.4889286723129618, 0.0000000000000002
+, 0.1942962641619352, 0.0000000000000001, 0.3158440934858354, 0.0
+000000000000001, 0.2496900204706467, 0.0000000000000001, 0.2959512
+116085838, 0.0000000000000000, 0.1552508901502876, 0.0000000000000
+001, 0.3427579648786467, 0.0000000000000002, 0.0541435546335990,
+0.0000000000000000, 0.2876031326588631, 0.0000000000000001, 0.2927
+760031123848 ],
[ 0.0000181563100889, 0.0000000000000002, 0.0005488425504935, 
+0.0000000000000000, 0.0060688055327925, 0.0000000000000002, 0.0372
+631786179999, 0.0000000000000002, 0.1425926143202678, 0.0000000000
+000000, 0.3478710011662555, 0.0000000000000002, 0.5076314191518138
+, 0.0000000000000006, 0.3123574515561076, 0.0000000000000001, 0.2
+091284023204721, 0.0000000000000002, 0.3689178831666572, 0.0000000
+000000000, 0.1344053015043881, 0.0000000000000001, 0.3299092963194
+629, 0.0000000000000004, 0.2082278951025315, 0.0000000000000003, 
+0.2240489250735314, 0.0000000000000001, 0.3188632967650468, 0.0000
+000000000002 ],
[ 0.0000000000000000, 0.0000646482239263, 0.0000000000000000,
+0.0011385375166878, 0.0000000000000000, 0.0098338922987181, 0.0000
+000000000000, 0.0513935599025939, 0.0000000000000001, 0.1741934371
+091206, 0.0000000000000002, 0.3836429445985159, 0.0000000000000004
+, 0.5045604175490942, 0.0000000000000004, 0.2522551752686503, 0.0
+000000000000001, 0.2579441025360375, 0.0000000000000004, 0.3397702
+320139469, 0.0000000000000003, 0.1752821352580565, 0.0000000000000
+001, 0.3071408379256653, 0.0000000000000001, 0.2279977492913131,
+0.0000000000000002, 0.2051227648591225, 0.0000000000000004, 0.3168
+989287342151 ],
[ 0.0000019661699240, 0.0000000000000002, 0.0000719533317298,
+0.0000000000000001, 0.0009846334801556, 0.0000000000000002, 0.0077
+074332169402, 0.0000000000000004, 0.0392493847315142, 0.0000000000
+000002, 0.1365860267496776, 0.0000000000000000, 0.3260488848298075
+, 0.0000000000000001, 0.5075095625005354, 0.0000000000000000, 0.4
+189127624432365, 0.0000000000000003, 0.0412751650984030, 0.0000000
+000000002, 0.4061413387987063, 0.0000000000000002, 0.1232820067450
+982, 0.0000000000000001, 0.3525839817655427, 0.0000000000000002,
+0.1125603524749182, 0.0000000000000002, 0.3486308267332468, 0.0000
+000000000004 ],
[ 0.0000000000000000, 0.0000071522511841, 0.0000000000000000, 
+0.0001531738843782, 0.0000000000000000, 0.0016457996800569, 0.0000
+000000000000, 0.0110351632854802, 0.0000000000000000, 0.0502042359
+369793, 0.0000000000000001, 0.1598827125229744, 0.0000000000000002
+, 0.3541055995816644, 0.0000000000000003, 0.5131002961113091, 0.0
+000000000000004, 0.3827211429896542, 0.0000000000000001, 0.0868953
+534697583, 0.0000000000000002, 0.4060663544932935, 0.0000000000000
+002, 0.0906849061585278, 0.0000000000000001, 0.3565218256839784,
+0.0000000000000002, 0.0954798894669647, 0.0000000000000000, 0.3460
+718461864576 ],
[ 0.0000001292972631, 0.0000000000000001, 0.0000057407224638,
+0.0000000000000000, 0.0000969677496746, 0.0000000000000000, 0.0009
+576391156081, 0.0000000000000002, 0.0063336604552633, 0.0000000000
+000000, 0.0298194086989161, 0.0000000000000001, 0.1025898937974792
+, 0.0000000000000004, 0.2581306279126550, 0.0000000000000004, 0.4
+608777058952639, 0.0000000000000003, 0.5289025115880573, 0.0000000
+000000000, 0.2501979027713430, 0.0000000000000003, 0.2459229288491
+326, 0.0000000000000000, 0.3976890721949536, 0.0000000000000001, 
+0.0441899989974140, 0.0000000000000001, 0.3828591066315501, 0.0000
+000000000001 ],
[ 0.0000000000000000, 0.0000004836289753, 0.0000000000000000,
+0.0000125935691119, 0.0000000000000000, 0.0001674328686760, 0.0000
+000000000000, 0.0014204630187498, 0.0000000000000000, 0.0084233556
+059741, 0.0000000000000000, 0.0364726872642141, 0.0000000000000000
+, 0.1173353239240527, 0.0000000000000001, 0.2791623093846674, 0.0
+000000000000001, 0.4743110004921181, 0.0000000000000002, 0.5168194
+754141317, 0.0000000000000002, 0.2172212810095154, 0.0000000000000
+001, 0.2660675030932704, 0.0000000000000000, 0.3868615369372250, 
+0.0000000000000001, 0.0582093189941330, 0.0000000000000002, 0.3792
+058820414961 ],
[ 0.0000000043467313, 0.0000000000000002, 0.0000002354338414,
+0.0000000000000001, 0.0000049191205605, 0.0000000000000001, 0.0000
+611103982443, 0.0000000000000001, 0.0005191484660093, 0.0000000000
+000002, 0.0032251377193584, 0.0000000000000002, 0.0151840473609755
+, 0.0000000000000001, 0.0551258405915315, 0.0000000000000001, 0.1
+547931322781031, 0.0000000000000001, 0.3318050095086971, 0.0000000
+000000004, 0.5216362528850688, 0.0000000000000001, 0.5386522451902
+931, 0.0000000000000004, 0.2190409351878887, 0.0000000000000002, 
+0.2724800121518238, 0.0000000000000002, 0.4221250094454858, 0.0000
+000000000001 ],
[ 0.0000000000000000, 0.0000000167699136, 0.0000000000000000, 
+0.0000005331602679, 0.0000000000000000, 0.0000087760999773, 0.0000
+000000000000, 0.0000937543740401, 0.0000000000000000, 0.0007150017
+363018, 0.0000000000000000, 0.0040915831845336, 0.0000000000000000
+, 0.0180530719106190, 0.0000000000000000, 0.0621905012332527, 0.0
+000000000000000, 0.1672502296997740, 0.0000000000000001, 0.3457819
+876105404, 0.0000000000000002, 0.5268406404622072, 0.0000000000000
+002, 0.5273801655434678, 0.0000000000000002, 0.1998433947965230,
+0.0000000000000000, 0.2805004874860440, 0.0000000000000003, 0.4170
+641909080779 ],
[ 0.0000000000469012, 0.0000000000000001, 0.0000000031569888, 
+0.0000000000000001, 0.0000000829508766, 0.0000000000000000, 0.0000
+013139183774, 0.0000000000000001, 0.0000144662736142, 0.0000000000
+000001, 0.0001187954144757, 0.0000000000000001, 0.0007576595856397
+, 0.0000000000000002, 0.0038457059058298, 0.0000000000000001, 0.0
+157544770583618, 0.0000000000000002, 0.0524034914224598, 0.0000000
+000000002, 0.1412730634268711, 0.0000000000000000, 0.3051191205191
+895, 0.0000000000000000, 0.5135323458897670, 0.0000000000000002,
+0.6319248579965772, 0.0000000000000001, 0.4699894042726668, 0.0000
+000000000000 ],
[ 0.0000000000000000, 0.0000000001872384, 0.0000000000000000,
+0.0000000073974097, 0.0000000000000000, 0.0000001531161950, 0.0000
+000000000000, 0.0000020854332860, 0.0000000000000000, 0.0000206106
+783438, 0.0000000000000000, 0.0001558996179038, 0.0000000000000000
+, 0.0009318666076831, 0.0000000000000000, 0.0044889011562547, 0.0
+000000000000000, 0.0176207946824076, 0.0000000000000000, 0.0565953
+635430333, 0.0000000000000000, 0.1482728021860708, 0.0000000000000
+000, 0.3129386149972143, 0.0000000000000000, 0.5172629611893457, 
+0.0000000000000001, 0.6280792712999990, 0.0000000000000001, 0.4631
+473700646182 ],
] ); # matrix T de los vectores del doblepozo
my $vect = transpose( $vec ); # matrix con los vectores en la columna
# <nVm>
my $V = $vec x $X x $vect;
my $step = 0.1;
my $hbar = 1;
#print Vnm( 0,1 );
#print Vnm( 1,0 );
my @Em = (
3.49090786, 3.48609059, 0.94249804,
+ 0.72223315, 0.84486159,
1.98334101, 3.52168993, 5.19562490,
+ 7.00900149, 8.96758961,
11.0614444444542760, 13.1282317263938650, 15.4597663541535190,
+ 18.2923955058715430, 19.5751133106800170,
21.3254642902799140, 27.3060148356095880, 33.6153274807702470,
+ 48.2864949409823790, 58.1657460989797170,
83.0782032408711330, 97.3475523518670660, 137.5307666823436800,
+157.3985331875067300, 221.5161226004345300,
248.6116702653141900, 353.3740736965299900, 390.0680354316574500,
+577.9792266930371600, 628.4681868621187300,
);
my $w = $Em[ 1 ]  $Em[ 0 ];
sub Vnm {
my ( $x1, $x2 ) = @_;
return $V>range( [ $x1, $x2 ] );
}
my @c = ( 1, 0, 0, 0 ); #INITIAL CONDITIONS
#@c = F = ( F0, F1, F2, F3 )
sub F {
my ( $m, $t, @c ) = @_;
my $f = 0;
for my $n ( 0 .. 3 ) {
my $a = Vnm( $n, $m ) * exp( ( i ) * ( $Em[ $m ]  $Em[ $n ]
+) * ( $t / $hbar ) );
if ( $n != $m ) {
$f += $a * ( ( i ) * ( $eprim ) * cos( $w * $t ) * $c[ $
+n ] );
}
}
return $f;
}
sub RK4_c {
my ( $m, $t, @c ) = @_;
my $k1 = sub {
my ( $t, @c ) = @_;
return $step * F( $m, $t, @c );
};
my $k2 = sub {
my ( $t, @c ) = @_;
return $step * F( $m, ( $t + $step / 2 ), ( @c + ( $k1>( $t )
+ / 2 ) ) );
};
my $k3 = sub {
my ( $t, @c ) = @_;
return $step * F( $m, ( $t + $step / 2 ), ( @c + ( $k2>( $t )
+ / 2 ) ) );
};
my $k4 = sub {
my ( $t, @c ) = @_;
return $step * F( $m, ( $t + $step ), ( @c + $k3>( $t )
+ ) );
};
return $c[ $m ] + ( 1 / 6 ) * ( $k1>( $t, @c ) + 2 * ( $k2>( $t,
+ @c ) + $k3>( $t, @c ) ) + $k4>( $t, @c ) );
}
open my $FILE, '>', 'C0_V1W01p2.dat';
for ( my $t = 0; $t < 2500; $t += $step ) {
#print " @c $t \n";
#print abs( $c[ 0 ] ) * abs( $c[ 0 ] ), "\n";
my ( $x, $b, $s, $d ) = map abs( RK4_c( $_, $t, @c ) ) ** 2, 0 ..
+3;
print $FILE "$x $b $s $d $t\n";
}
close $FILE;
 [reply] [d/l] 
Re: RK4 for a vectorial function
by bliako (Monsignor) on Nov 26, 2021 at 08:33 UTC

for 0<t<25 I get beautiful oscillations in columns 3 and 4 of the output file
side remark: You don't need to declare subs k1_1 etc. as variables (e.g. my $k1_1;) (and you could also implement them outside the RK4_... subs, EDIT: did you mean to use them as closures? you don't seem to because you pass them the important params and step is "global" in their scope). EDIT: also you can stop calling subs using the &mysub() notation, instead mysub() is valid.
 [reply] [d/l] [select] 
Re: RK4 for a vectorial function
by etj (Pilgrim) on Dec 18, 2021 at 15:18 UTC

Sorry to have only just noticed this; I see you're using PDL::Complex, which is now officially deprecated. If you delete that line, then with a recent (2.047+) PDL "i" will return a "cdouble" (a "native complex" type) with the obvious value and all the type stuff will do the right thing.  [reply] 

