# Functional Programming for Non-linear Equations

GROUPS:
 Pat Mac 1 Vote Hello All,I am on a journey to become adept at developing functional programs with Mathematica. I have a short prototypecode that I am using to learn how to convert procedural programs into functional programs and am wondering if anyonecould provide any more advise? Specifically I would like to learn how to write this code in a much more memory friendly manner.And also learn a better way to run the main loop of the Markov Chain Monte Carlo simulation that minimizes the distance between the experimental radioactive decay data and my model. Does anyone have any suggestions on this or how to develop into better functional programmer.Thanks,Pat  Clear["Global`*"]    texp = {0.2, 2.2, 4.0, 5.0, 6.0, 8.0, 11.0, 15.0, 18.0, 26.0, 33.0,      39.0, 45.0};  dexp = {35., 25., 22.1, 17.9, 16.8, 13.7, 12.4, 7.5, 4.9, 4.0, 2.4,      1.4, 1.1};  dt = Transpose[{texp, dexp}];  dataPlot = ListPlot[     dt,    PlotStyle -> {Red, AbsolutePointSize[6]},    Frame -> {True, True, False, False}, Axes -> False,    FrameLabel -> {"Time", "Species A"},    LabelStyle -> {FontSize -> 20, FormatType -> Bold},    RotateLabel -> False,    ImageSize -> {550}];  runlength = 10000; sets = runlength/10.; testing = 10.0; Mean[dexp]; StandardDeviation[dexp]; stim = 0.01;  chi2[ka_?NumberQ, kb_?NumberQ] := Block[{sol, A, B, c}, sol =     NDSolve[       { A'[t] == -ka*A[t],        B'[t] == ka*A[t] - kb*B[t],        c'[t] == kb*B[t],        A[0] == 35.,        B[0] == 0.,        c[0] == 0.},       {A, B, c},       {t, 0., 100.}][[1]];    Apply[Plus, (dexp - (A[t] /. sol /. t -> texp))^2]];  DateString[]; RandomSeed[AbsoluteTime[%]]; (* Set Random Seed based on current time *)  Clear[deltaE] Clear[freeE] Clear[ka] Clear[kb] Clear[u, v] bestParamSets = {}; thrownsets = {}; (* These arrays are for storage *)  ka = ConstantArray[0, {runlength}]; kb = ConstantArray[0, {runlength}]; deltaE = ConstantArray[0, {runlength}]; freeE = ConstantArray[0, {runlength}]; v = ConstantArray[0, {runlength}]; u = ConstantArray[0, {runlength}]; decayscore = ConstantArray[0, {runlength}]; (**************************************************) (******** Starting the Metropolis Algorithm ********) (**************************************************)  range = 10.0;(* Pick different initial guesses for param. sets *) weight = 1./10; ka[[0]] = RandomReal[{0.0, 0.2}]; kb[[0]] = RandomReal[{0.0, 0.2}]; least = 0.0004; max = 1.000; decayscore[[0]] = chi2[ka[[0]], kb[[0]]]; beta = 1/.01;   (tbl = Reap[      Table[(*Open Main For Loop*)   ka[[i]] = ka[[i - 1]] + weight*ka[[i - 1]]*RandomReal[{-1.0, 1.0}];   kb[[i]] = kb[[i - 1]] + weight*kb[[i - 1]]*RandomReal[{-1.0, 1.0}];       (*************************************************)       (*This is where the evaluation of the parameter sets begins*)       (*************************************************)      If[      least < ka[[i]] && ka[[i]] < max &&      least < kb[[i]] && kb[[i]] < max,      decayscore[[i]] = chi2[ka[[i]], kb[[i]]];(*This is the function evaluated for each parameter set*)      deltaE[[i]] = decayscore[[i]] - decayscore[[i - 1]];(*The change in energy for each step*)      freeE[[i]] = Exp[-beta*deltaE[[i]]];(*This is a free energy associated with each step*)        v[[i]] = Min[1, freeE[[i]]];        u[[i]] = RandomReal[{0, 1}];        (**************************************************)        (******** Metropolis Selection Algorithm ********)        (**************************************************)      If[u[[i]] < v[[i]],(*If u {data}]; ) // AbsoluteTiming  (*Close Main For Loop How long does this loop take to run?*)
5 years ago
4 Replies
 Is this question ill posed? I've noticed that no one has responded.