2
$\begingroup$

I have the following code generated in Mathematica 7. When I DSolve it in Mathematica 7, the imaginary parts of the solution are really low-valued. As a result, I can easily chop them and plot the results. However, the same operation does not give low values for the imaginary parts in Mathematica 10.4 and the results cannot be plotted.

Here is the code.

Ef = 220000; bf = 1; tf = 0.167; Af = bf*tf;

EA = Ef*Af; EI = 100; kt = 10; kr = 0.1; r = 25; GA = 
 100000*Af; \[Chi] = 1.2;

ut[s_];
ur[s_];

pt[s_] = kt*ut[s];
pr[s_] = kr*ur[s];

equ1 = n'[s] - T[s]/r + pt[s];
equ2 = T'[s] + n[s]/r + pr[s];
equ3 = M'[s] - T[s];
eps[s_] = ut'[s] - ur[s]/r;
\[Gamma]s[s_] = ur'[s] + ut[s]/r + \[Phi][s];
k[s_] = \[Phi]'[s];
equ4 = eps[s] - n[s]/EA;
equ5 = k[s] + M[s]/EI;
equ6 = \[Gamma]s[s] - (\[Chi]*T[s])/GA;

an = DSolve[{equ1 == 0, equ2 == 0, equ3 == 0, equ4 == 0, equ5 == 0, 
    equ6 == 0,
    n[0] - F == 0, n[\[Pi]/2*r] - F == 0,
    ut[0] == 0, ur[0] == 0, ut[\[Pi]/2*r] == 0, 
    ur[\[Pi]/2*r] == 0}, {n[s], T[s], M[s], ut[s], ur[s], \[Phi][s]}, 
   s];

n[s_] = Chop[n[s] /. an[[1, 2]]];

Plot[n[s], {s, 0, \[Pi]/2*r}]

The results of different versions of Mathematica for the solution of the differential equation is:

Mathematica 7:

{{n[s] -> 
   E^(-0.345376 s) ((1.19981 - 1.91926*10^-16 I) E^(0.128248 s)
        F + (0.000237717 + 9.04032*10^-19 I) E^(0.562504 s)
        F + (0.00798872 + 2.12516*10^-17 I) E^(0.473623 s)
        F Cos[0.204351 s] - (0.208036 - 1.6977*10^-16 I) E^(
       0.217128 s)
        F Cos[0.204351 s] + (6.72924*10^-6 - 2.64133*10^-17 I) E^(
       0.473623 s)
        F Sin[0.204351 s] + (1.21171 - 1.21342*10^-16 I) E^(
       0.217128 s) F Sin[0.204351 s]), 
  T[s] -> E^(-0.345376 s) ((0.153989 - 2.46326*10^-17 I) E^(
       0.128248 s) F - (0.0000305097 + 1.16027*10^-19 I) E^(
       0.562504 s) F - (0.000441035 + 9.1325*10^-18 I) E^(0.473623 s)
        F Cos[0.204351 s] + (0.217962 - 2.29161*10^-17 I) E^(
       0.217128 s)
        F Cos[0.204351 s] - (0.00151332 + 1.9483*10^-17 I) E^(
       0.473623 s)
        F Sin[0.204351 s] + (0.106487 - 2.85166*10^-17 I) E^(
       0.217128 s) F Sin[0.204351 s]), 
  M[s] -> E^(-0.345376 s) ((-0.709208 + 1.13448*10^-16 I) E^(
       0.128248 s) F - (0.000140515 + 5.34374*10^-19 I) E^(0.562504 s)
        F + (0.0043412 + 2.14946*10^-17 I) E^(0.473623 s)
        F Cos[0.204351 s] - (0.85409 - 2.14849*10^-16 I) E^(
       0.217128 s)
        F Cos[0.204351 s] - (0.00488269 + 7.14507*10^-17 I) E^(
       0.473623 s)
        F Sin[0.204351 s] + (0.53059 + 1.99297*10^-17 I) E^(
       0.217128 s) F Sin[0.204351 s]), 
  ut[s] -> E^(-0.345376 s) ((0.0266672 - 4.26578*10^-18 I) E^(
       0.128248 s) F - (5.28354*10^-6 + 2.00932*10^-20 I) E^(
       0.562504 s) F - (0.000104355 + 1.08889*10^-18 I) E^(0.473623 s)
        F Cos[0.204351 s] - (0.0265575 - 5.37477*10^-18 I) E^(
       0.217128 s)
        F Cos[0.204351 s] + (0.000157111 + 4.80885*10^-19 I) E^(
       0.473623 s)
        F Sin[0.204351 s] + (0.0117146 - 3.91483*10^-19 I) E^(
       0.217128 s) F Sin[0.204351 s]), \[Phi][s] -> 
   E^(-0.345376 s) ((-0.0326632 + 5.22492*10^-18 I) E^(0.128248 s)
        F + (6.47152*10^-6 + 2.4611*10^-20 I) E^(0.562504 s)
        F - (0.00026707 + 2.52378*10^-18 I) E^(0.473623 s)
        F Cos[0.204351 s] - (0.000190439 - 3.34787*10^-18 I) E^(
       0.217128 s)
        F Cos[0.204351 s] - (0.0000448289 - 1.24984*10^-18 I) E^(
       0.473623 s)
        F Sin[0.204351 s] + (0.0416757 - 5.1102*10^-18 I) E^(
       0.217128 s) F Sin[0.204351 s]), 
  ur[s] -> E^(-0.345376 s) ((-0.145571 + 2.32861*10^-17 I) E^(
       0.128248 s) F - (0.0000288419 + 1.09685*10^-19 I) E^(
       0.562504 s) F + (0.000462626 - 6.46096*10^-18 I) E^(0.473623 s)
        F Cos[0.204351 s] + (0.145137 - 1.67154*10^-17 I) E^(
       0.217128 s)
        F Cos[0.204351 s] + (0.00103685 + 2.53644*10^-18 I) E^(
       0.473623 s)
        F Sin[0.204351 s] + (0.097293 - 2.50601*10^-17 I) E^(
       0.217128 s) F Sin[0.204351 s])}}

Mathematica 10.4:

{{M[s] -> (3.52596 + 
      3.28357 I) E^(-0.345376 s) ((0. - 1. I) E^(0.128248 s)
        F + (0.00444174 + 0.0871541 I) E^(0.562504 s)
        F - (0.392686 + 0.29177 I) E^(0.473623 s)
        F Cos[0.204351 s] + (0.252446 - 1.54495 I) E^(0.217128 s)
        F Cos[0.204351 s] + (1.60583 + 0.145085 I) E^(0.473623 s)
        F Sin[0.204351 s] + (1.90461 - 0.0680882 I) E^(0.217128 s)
        F Sin[0.204351 s]), 
  n[s] -> (5.96507 + 
      5.55501 I) E^(-0.345376 s) ((0. + 1. I) E^(0.128248 s)
        F - (0.00444174 + 0.0871541 I) E^(0.562504 s)
        F - (1.05642 + 0.218411 I) E^(0.473623 s)
        F Cos[0.204351 s] + (1.15065 - 0.778044 I) E^(0.217128 s)
        F Cos[0.204351 s] + (0.558183 - 0.0881049 I) E^(0.473623 s)
        F Sin[0.204351 s] + (0.777341 + 0.802445 I) E^(0.217128 s)
        F Sin[0.204351 s]), 
  T[s] -> (0.765583 + 
      0.712954 I) E^(-0.345376 s) ((0. + 1. I) E^(0.128248 s)
        F + (0.00444174 + 0.0871541 I) E^(0.562504 s)
        F + (1.27939 - 0.0357876 I) E^(0.473623 s)
        F Cos[0.204351 s] + (1.64343 + 0.848452 I) E^(0.217128 s)
        F Cos[0.204351 s] + (1.31807 + 0.360296 I) E^(0.473623 s)
        F Sin[0.204351 s] - (1.36256 - 1.49426 I) E^(0.217128 s)
        F Sin[0.204351 s]), 
  ur[s] -> (0.723733 + 
      0.673981 I) E^(-0.345376 s) ((0. - 1. I) E^(0.128248 s)
        F + (0.00444174 + 0.0871541 I) E^(0.562504 s)
        F - (1.10206 + 0.0102304 I) E^(0.473623 s)
        F Cos[0.204351 s] + (1.38147 + 0.486021 I) E^(0.217128 s)
        F Cos[0.204351 s] - (0.862745 + 0.275686 I) E^(0.473623 s)
        F Sin[0.204351 s] - (0.858699 - 1.21572 I) E^(0.217128 s)
        F Sin[0.204351 s]), 
  ut[s] -> (0.132581 + 
      0.123467 I) E^(-0.345376 s) ((0. + 1. I) E^(0.128248 s)
        F + (0.00444174 + 0.0871541 I) E^(0.562504 s)
        F + (0.125917 + 0.206204 I) E^(0.473623 s)
        F Cos[0.204351 s] - (0.0128045 + 1.16713 I) E^(0.217128 s)
        F Cos[0.204351 s] - (1.26293 + 0.141651 I) E^(0.473623 s)
        F Sin[0.204351 s] + (1.47499 - 0.217812 I) E^(0.217128 s)
        F Sin[0.204351 s]), \[Phi][
    s] -> (0.162391 + 
      0.151227 I) E^(-0.345376 s) ((0. - 1. I) E^(0.128248 s)
        F - (0.00444174 + 0.0871541 I) E^(0.562504 s)
        F + (1.41196 + 0.250179 I) E^(0.473623 s)
        F Cos[0.204351 s] + (1.57263 - 0.791006 I) E^(0.217128 s)
        F Cos[0.204351 s] - (0.468887 - 0.153003 I) E^(0.473623 s)
        F Sin[0.204351 s] + (0.71873 + 1.14512 I) E^(0.217128 s)
        F Sin[0.204351 s])}}

As it is obvious the imaginary parts for the result of Mathematica 7 are negligible. While in the case of Mathematica 10.4 Chop does not do anything and the imaginary parts are significant.

$\endgroup$
  • $\begingroup$ Can you provide the different outputs? $\endgroup$ – Sektor May 22 '17 at 15:23
  • $\begingroup$ You did not include a value for F so your code cannot be evaluated. Rationalize all of the inexact numbers prior to DSolve, Using F=1; the Chop is not necessary and the plot evaluates successfully with version 10.4.1. $\endgroup$ – Bob Hanlon May 22 '17 at 15:25
  • $\begingroup$ The post is updated. I am intended to obtain the same solution as in Mathematica 7 but I can't. $\endgroup$ – KratosMath May 22 '17 at 15:30
  • $\begingroup$ @BobHanlon Thanks a lot for your answer. Yes, that did solve my problem. $\endgroup$ – KratosMath May 22 '17 at 16:55
2
$\begingroup$
$Version

(*  "10.4.1 for Mac OS X x86 (64-bit) (April 11, 2016)"  *)

Clear[n]
Ef = 220000; bf = 1; tf = 0.167; Af = bf*tf;

EA = Ef*Af; EI = 100; kt = 10; kr = 0.1; r = 25; GA = 100000*Af; χ = 1.2;

ut[s_];
ur[s_];

pt[s_] = kt*ut[s];
pr[s_] = kr*ur[s];

equ1 = n'[s] - T[s]/r + pt[s];
equ2 = T'[s] + n[s]/r + pr[s];
equ3 = M'[s] - T[s];
eps[s_] = ut'[s] - ur[s]/r;
γs[s_] = ur'[s] + ut[s]/r + ϕ[s];
k[s_] = ϕ'[s];
equ4 = eps[s] - n[s]/EA;
equ5 = k[s] + M[s]/EI;
equ6 = γs[s] - (χ*T[s])/GA;

Imaginary artifacts and need for Chop are due to precision issues. Use exact numbers (Rationalize) to avoid loss of precision.

Clear[n];
an = DSolve[{equ1 == 0, equ2 == 0, equ3 == 0, equ4 == 0, equ5 == 0, equ6 == 0,
       n[0] - F == 0, n[π/2*r] - F == 0, ut[0] == 0, ur[0] == 0, 
      ut[π/2*r] == 0, ur[π/2*r] == 0} // Rationalize, {n[s], T[s], 
     M[s], ut[s], ur[s], ϕ[s]}, s] // Simplify;

n[s_] = n[s] /. an[[1, 2]];

n[1.]

(*  1.01194 F  *)

F needs to have a value to Plot the value of n[s]

Plot[n[s] /. F -> 1, {s, 0, π/2*r}]

enter image description here

$\endgroup$
  • $\begingroup$ great solution!!! $\endgroup$ – KratosMath May 22 '17 at 17:14

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.