I am having an issue with the functions that I have defined using the solutions given by NDSolve. The easiest way to explain my issue is through the code itself, so here it is:

boundary1 = 0.3; 
boundary2 = 1; 

h[r_] = 10^(-2); 
\[Mu]1 = 8*10^4; 
\[Mu]2 = 8*10^4; 
p = 1300; 
nsolution1[(w_)?NumericQ] := NDSolve[{lt[r] + r*Derivative[1][lt][r] == lr[r]*Cos[\[Beta][r]], p*r*lt[r]*lr[r]*Cos[\[Beta][r]] - \[Mu]1*h[r]*Sin[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) - \[Mu]1*h[r]*r*Derivative[1][\[Beta]][r]*Cos[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) - 
        \[Mu]1*h[r]*r*Sin[\[Beta][r]]*(Derivative[1][lr][r] + (3*Derivative[1][lr][r])/(lr[r]^4*lt[r]^2) + (2*Derivative[1][lt][r])/(lr[r]^3*lt[r]^3)) == 0, (-p)*r*lr[r]*lt[r]*Sin[\[Beta][r]] - \[Mu]1*h[r]*Cos[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) + 
        \[Mu]1*h[r]*r*Derivative[1][\[Beta]][r]*Sin[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) - \[Mu]1*h[r]*r*Cos[\[Beta][r]]*(Derivative[1][lr][r] + (3*Derivative[1][lr][r])/(lr[r]^4*lt[r]^2) + (2*Derivative[1][lt][r])/(lr[r]^3*lt[r]^3)) + \[Mu]1*h[r]*(lt[r] - 1/(lr[r]^2*lt[r]^3)) == 0, lt[0.001] == w, 
      lr[0.001] == w, \[Beta][0.001] == 0.001}, {lr, lt, \[Beta]}, {r, 0.001, 1}][[1]]; 
radialstretch[(w_)?NumericQ, (r_)?NumericQ] := lr[r] /. nsolution1[w]; 
thetastretch[(w_)?NumericQ, (r_)?NumericQ] := lt[r] /. nsolution1[w]
angle[(w_)?NumericQ, (r_)?NumericQ] := \[Beta][r] /. nsolution1[w]
radialstretchderivative[(w_)?NumericQ, (r_)?NumericQ] := D[lr[u] /. nsolution1[w], u] /. u -> r
thetastretchderivative[(w_)?NumericQ, (r_)?NumericQ] := D[lt[u] /. nsolution1[w], u] /. u -> r
R[(w_)?NumericQ, (r_)?NumericQ] := r*thetastretch[w, r]
Rderivative[(w_)?NumericQ, (r_)?NumericQ] := thetastretch[w, r] + r*thetastretchderivative[w, r]
f1[(w_)?NumericQ, (r_)?NumericQ] := -(radialstretch[w, r]^2 - Rderivative[w, r]^2)^(1/2)
Z[(w_)?NumericQ, (r_)?NumericQ] := NIntegrate[Re[f1[w, t]], {t, 0.001, r}]

When I try to evaluate any of the functions I've defined at specific values of r, I do get an answer for them, however, when I try plotting the functions, that does not work. I understand it has to do with the way Mathematica is evaluating everything, however, I do not understand it well enough to know where exactly the problem is. Here is an example:

Plot[Evaluate[lr[r] /. nsolution1[1.2]], {r, 0.001, 1}]

This gives me a plot that is correct, but if I do the following:

Plot[radialstretch[1.2, r], {r, 0.001, 1}]

which should give me the same exact plot, instead I get the error: "NDSolve::dsvar: 0.0010204081428571428` cannot be used as a variable". together with a bunch of other ones.

Also, I thought either ParametricNDSolve or ParametricNDSolveValue might be better options in my case since I have w as a variable, however, that is somehow giving me a different answer for the functions when compared to the answers using NDSolve. I might have to make a separate question for that because I am extremely confused as to why I would get a different answer.


1 Answer 1


Look at the following:

x = 1.; NDSolve[{f'[x] == f[x], f[0] == 1}, f, {x, 0, 1}]

(* NDSolve::dsvar: 1.` cannot be used as a variable. *)


NDSolve[{f'[x] == f[x], f[0] == 1}, f, {x, 0, 1}];

Shows that x is defined as a global variable. However, dummy variables should be localized. I think this is a bug. A dirty fix is to use another variable in the plot:

Plot[radialstretch[1.2, rr], {rr, 0.001, 1}]

enter image description here


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.