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.

Improve this question
$\endgroup$
4
  • $\begingroup$ Can you provide the different outputs? $\endgroup$
    – Sektor
    Commented May 22, 2017 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
    Commented May 22, 2017 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
    Commented May 22, 2017 at 15:30
  • $\begingroup$ @BobHanlon Thanks a lot for your answer. Yes, that did solve my problem. $\endgroup$
    – KratosMath
    Commented May 22, 2017 at 16:55

1 Answer 1

Reset to default
4
$\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

Improve this answer
$\endgroup$
1
  • $\begingroup$ great solution!!! $\endgroup$
    – KratosMath
    Commented May 22, 2017 at 17:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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