How can I rewrite an equation written in FORTRAN file by Matlab?

Please, I have an equation differentiated in Matlab and transferred to a FORTRAN file. How can I get back the original equation (before differentiation) as the result written in FORTRAN file looks difficult to be arranged. Below is an example of the differentiated equation written in FORTRAN.

      t3 = EX*6.412498652425423D23
      t2 = t3-1.584301877890218D19
      t4 = t2**2
      t5 = t4+1.097400831008096D40
      t6 = EX*8.224027793471572D47
      t7 = t6-2.031866731401218D43
      t8 = 1.0D0/sqrt(t5)
      t9 = t7*t8*(1.0D0/2.0D0)
      t10 = t9+6.412498652425423D23
      t11 = sqrt(t5)
      t12 = t3+t11-1.584301877890218D19
      A0(1,1) = t10*1.0D0/t12**(2.0D0/3.0D0)*(1.0D0/3.0D0)+t10*1.0D0/t12
     +**(4.0D0/3.0D0)*7.407423474671637D12

[Walter Roberson]

With the symbolic toolbox it is not difficult to put that together into a single expression for A0 in terms of EX, getting

2.666666666666667*(0.7479709705395701e36*EX-0.1847972026918085e32+4777683803756105.*(0.2450951276264063e41*EX^2-0.1211086947561036e37*EX+0.6690627070730690e33)^(1/2))*(0.1025999784388068e26*EX+65536.*(0.2450951276264063e41*EX^2-0.1211086947561036e37*EX+0.6690627070730690e33)^(1/2)-0.2534883004624349e21+5688901228547817.*(0.1565551428814800e21*EX+(0.2450951276264063e41*EX^2-0.1211086947561036e37*EX+0.6690627070730690e33)^(1/2)-3867924506567915.)^(1/3))/((0.2450951276264063e41*EX^2-0.1211086947561036e37*EX+0.6690627070730690e33)^(1/2)*(0.1565551428814800e21*EX+(0.2450951276264063e41*EX^2-0.1211086947561036e37*EX+0.6690627070730690e33)^(1/2)-3867924506567915.)^(5/3))

However, there is no obvious integral for that.

[ND]

Thank you Walter,

Please how can I put the differentiated equation (A0) in a single expression?

Doing the integration would not be a problem, it can be done manually.

Answer 1

[Walter Roberson]

Rewrite it slightly. Change the D to E. Change the to ^ .

   syms EX

and then execute the code.

syms EX
t3 = EX*6.412498652425423E23
      t2 = t3-1.584301877890218E19
      t4 = t2^2
      t5 = t4+1.097400831008096E40
      t6 = EX*8.224027793471572E47
      t7 = t6-2.031866731401218E43
      t8 = 1.0E0/sqrt(t5)
      t9 = t7*t8*(1.0E0/2.0E0)
      t10 = t9+6.412498652425423E23
      t11 = sqrt(t5)
      t12 = t3+t11-1.584301877890218E19
      A0(1,1) = t10*1.0E0/t12^(2.0E0/3.0E0)*(1.0E0/3.0E0)+t10*1.0E0/t12^(4.0E0/3.0E0)*7.407423474671637E12

The result will be an expression stored in A0. All of the numeric values will have been converted to fractions, so it will start like

((822402779347157244357116847336510700427459690496*EX - 20318667314012179466854033896368807794966528)/(8192*((641249865242542279229440*EX

For numeric integration, you could either use vpaintegral on A0 or you could use matlabFunction(A0) and use integral() with the resulting function handle. Either way you would need numeric upper and lower bounds.

[Walter Roberson]

Oh yes, it was necessary to detect the continuation marker on the last line and merge with the previous line.

[ND]

I am nearly there, but how can I do the integral? could you please clarify this more?

Many thanks

[Walter Roberson]

For lower bound 2, upper bound 11:

vpaintegral(A0, 2, 11)

Or

F = matlabFunction(A0) ;
integral(F, 2, 11)

[ND]

Sorry I got number with this. How can I get the original equation after the integration?

[Walter Roberson]

What was the original equation that was differentiated to give the Fortran code? With that equation perhaps I can work backwards from the expression to see how you would have to integrate. But at the moment my belief is that the Fortran code is an approximation of the derivative rather than the whole derivative, such as if a Taylor expansion had been done.

[Walter Roberson]

I told you yesterday that the equation does not have any obvious closed form integral and you said that would not be a problem, so you should go ahead with whatever technique you had in mind at that time.

[Walter Roberson]

Numeric integration techniques have difficulty with negative EX starting at approximately -45/1000 (integrating to 0). For example by -51/1000 then using 50 digits is considered to be not enough to integrate to 0 without unacceptable loss of precision.

The equation is not quite symmetric around 0.

The equation is fairly steep near 0, but is not infinite: it peaks at about 1.904514780*10^10.