0
$\begingroup$

I am trying to solve a system $dx/dt =f (x, t)$ using the RK4 in Fortran 90 that works for any 2D Hamiltonian systems. I would like to generate using Mathematica a Fortran file .f90 that contains the terms (in FORTRAN SYNTAX) of the function $f(x, t)$ that corresponds to the right hand side of the RK4 solver. In the code it corresponds to the term deriv.

ndof = 2; ndim = 2*ndof;
Clear[h];
h = (1/2) (x[3]^2 + x[4]^2) + aa*x[1]^2/2 - bb*Sqrt[x[1]^2]^3/3 + 
 aa*x[2]^2/2 - bb*Sqrt[x[2]^2]^3/3 + aa*(x[1] - x[2])^2/2 - 
 bb*Sqrt[(x[1] - x[2])^2]^3/3;
derivq = Table[D[ h, x[i + ndof]]   , {i, 1, ndof}]
derivp = Table[-D[h, x[i]] , {i, 1, ndof}]
deriv  = Join[derivq, derivp]
filename = "1d_graphene_rk.f90"
$Output = filename; 
SetOptions[$Output]
Do[Print["dxdt(" <> ToString[i] <> ")=" <> 
    ToString[FortranForm[FullSimplify[deriv[[i]]]]]];, {i, 1, ndim}]
$Output = OutputStream["stdout", 1];
SetOptions[$Output]

it produces the following errors when trying to save deriv

Print::noopen: Cannot open 1d_graphene_rk.f90.
Print::noopen: Cannot open 1d_graphene_rk.f90.
Print::noopen: Cannot open 1d_graphene_rk.f90.
General::stop: Further output of Print::noopen will be suppressed 
 during this calculation.

How can I fix this code for it to work?

$\endgroup$
  • 1
    $\begingroup$ Try $Output = OpenWrite["1d_graphene_rk.f90"]; $\endgroup$ – Pillsy Apr 24 '18 at 12:07
2
$\begingroup$
Export["test.f90", 
 Table["dxdt(" <> ToString[i] <> ")=" <> 
   ToString[FortranForm[FullSimplify[deriv[[i]]]]], {i, 
   ndim}], "Text"]

FilePrint["test.f90"]

enter image description here

aside you may want to add Assumptions -> Element[{x[1], x[2]}, Reals] to Simplify ( that will turn Sqrt[x[1]^2] to Abs[x[1]] )

note if you really want to do "low level" file i/o you should work with Write , WriteString etc instead of messing with Print / $Output

$\endgroup$
  • $\begingroup$ Running the command as it is on the post, I am getting the following error OpenWrite::noopen: Cannot open /test.f90. '$Failed`. What is happening on my notebook :-( :-( ? $\endgroup$ – many Apr 25 '18 at 15:32
  • 1
    $\begingroup$ oh you may be running in a directory where you do not have write permission. Do use SetDirectory to set to somplace you can write or give a full path name to the file. $\endgroup$ – george2079 Apr 26 '18 at 0:32
  • $\begingroup$ Great, many thanks. It works fine now. $\endgroup$ – many Apr 26 '18 at 10:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.