2

I wanted to write this code in such a way that I run the whole thing for each $n$, and not write a1 a2 a3 or eq1 eq2 eq3 every time.

Clear[y, x, M, V, th, EI, PE]
y = a4*x^4 + a3*x^3 + a2*x^2;
th = D[y, x];
M = EI*D[y, {x, 2}];
V = EI*D[y, {x, 3}];
PE = EI/2*Integrate[D[y, {x, 2}]^2, {x, 0, L}] + P*(y /. x -> L);
Eq2 = D[PE, a2];
Eq1 = D[PE, a3];
Eq3 = D[PE, a4]
Sol = Solve[{Eq1 == 0, Eq2 == 0, Eq3 == 0}, {a3, a2, a4}];
y = y /. Sol[[1]];
th /. Sol[[1]];
FullSimplify[M /. Sol[[1]]]
FullSimplify[V /. Sol[[1]]]

Here's my attempt but it doesn't work:

Clear[a, y, x, M, V, th, EI, PE]
n = 4
y = Sum[{y^i}*a[i], {i, 0, n}]
th = D[y, x];
M = EI*D[y, {x, 2}];
V = EI*D[y, {x, 3}];
PE = EI/2*Integrate[D[y, {x, 2}]^2, {x, 0, L}] + P*(y /. x -> L)
Eq = Table [{D[PE, a[i]]}, {i, 2, n}]
Sol = Solve[{Eq[i] == 0}, {a[i]}];
y = y /. Sol[[1]];
th /. Sol[[1]];
FullSimplify[M /. Sol[[1]]]
FullSimplify[V /. Sol[[1]]]
flinty
  • 25,147
  • 2
  • 20
  • 86
Ali AlCapone
  • 141
  • 6

1 Answers1

3

Be careful not to define things recursively like y. Use Array[a, n] to produce the list {a[1],a[2],a[3],a[4]}. You also do not need some of the lists {...} like in your Table and Solve, and you can use ConstantArray to make your zero vector for the equations:

Remove["Global`*"]
n = 4;
y = Sum[x^i*a[i], {i, 0, n}]
th = D[y, x];
M = EI*D[y, {x, 2}];
V = EI*D[y, {x, 3}];
PE = EI/2*Integrate[D[y, {x, 2}]^2, {x, 0, L}] + P*(y /. x -> L)
Eq = Table[D[PE, a[i]], {i, 2, n}]
sol = Solve[Eq == ConstantArray[0, n - 1], Array[a, n]]
y /. sol
th /. sol
Eq /. sol
FullSimplify[M /. sol[[1]]]
FullSimplify[V /. sol[[1]]]

(* {{a[2] -> -((L P)/(2 EI)), a[3] -> P/(6 EI), a[4] -> 0}} *)
flinty
  • 25,147
  • 2
  • 20
  • 86
  • Thanx, a lot I'm really new in this do you know where can I learn these simple codes more? – Ali AlCapone Aug 15 '20 at 23:39
  • 1
    @AliAlCapone no problem. Please upvote/accept if you find it useful. You may be interested in this since you're new to Mathematica: https://mathematica.stackexchange.com/questions/18393/what-are-the-most-common-pitfalls-awaiting-new-users – flinty Aug 15 '20 at 23:41
  • I did Thanx :-) – Ali AlCapone Aug 15 '20 at 23:43
  • 1
    +1 You can simplify the Solve statement to sol = Solve[Eq == 0, Array[a, n]]. Solve will treat the 0 as if it were a vector. Also, the warning message can be avoided with sol=Solve[Eq==0, Rest@Array[a, n]] – Bob Hanlon Aug 15 '20 at 23:53
  • is it possible in Mathematica to write an Array for n a plot the result for different n, so that I can compare every Plot with the others to compare? it is easy to do in python or Matlab but I'm a little childe in Mathematica – Ali AlCapone Aug 16 '20 at 02:57
  • 1
    @AliAlCapone which variables do you want to plot for different n? – flinty Aug 16 '20 at 12:18
  • @flinty I want to plot M, V and Y – Ali AlCapone Aug 16 '20 at 15:08
  • 1
    @AliAlCapone Do you want to plot M vs. x, V vs. x, and Y vs. x for various n ? In which case you need to provide values for EI , L, and P. Also the solution applied to Y is has floating a[0] and a[1] because sol only provides values for some of the a[i], so it doesn't make sense to plot -((L P x^2)/(2 EI)) + (P x^3)/(6 EI) + a[0] + x a[1] – flinty Aug 16 '20 at 15:11
  • @flinty M vs. x, V vs. x, and Y vs. x for various n. and i change it a little because a[0] and a[1] are 0 Remove["Global`*"] n = 4; y = Sum[x^i*a[i], {i, 2, n}]; th = D[y, x]; M = EI*D[y, {x, 2}]; V = EI*D[y, {x, 3}]; PE = EI/2*Integrate[D[y, {x, 2}]^2, {x, 0, L}] + P*(y /. x -> L); Eq = Table[D[PE, a[i]], {i, 1, n}]; sol = Solve[Eq == 0, Array[a, n]] y /. sol th /. sol Eq /. sol; FullSimplify[M /. sol[[1]]] FullSimplify[V /. sol[[1]]] EI = 100; P = 1; L = 10; Plot[{V, -100, 100}, {x, 0, L}, AxesLabel -> {"x (m)", "M (N.m.)"}, PlotLegends -> {Style[" M", Bold, 12]}] – Ali AlCapone Aug 16 '20 at 16:11