125 Reputation

6 Badges

6 years, 122 days

MaplePrimes Activity

These are replies submitted by maple2015

@Carl Love 

Sorry about incomplete code (h:=1)

you can get result for constant E (say E=200e9) or variable E and little amount of L, however for variable E and large amunt of L it does not work. 

I think one of the reasons is piecewise function. Is it possible to use a simple equivalent function instead of it?


T := time(): 
M := 20: 
Digits := 30: 
L := 500: 
h := 1: 
R := (1/2)*h: 
nu := 0.3: 
E := 0.200e12: 
X := (int(E*(z+(1/2)*R), [z = y-(1/2)*R .. 0, y = -(1/2)*R .. (1/2)*R]))/(int(E, [z = y-(1/2)*R .. 0, y = -(1/2)*R .. (1/2)*R])): 
beta := Pi^2/L^2: R := (1/2)*h: G := E/(2*(1+nu)): 
phi := add(b[n]*y^n, n = 0 .. M): 
Eq := diff(phi, y$2)+(diff(E, y))*(diff(phi, y))/E+((diff(E, y$2))/E-((diff(E, y))/E)^2)*phi-2*beta*(1+nu)*(phi-1): 
st := [seq(coeftayl(Eq, y = 0, j), j = 0 .. M-2)]: 
for k to M-1 do 
b[k+1] := solve(st[k], b[k+1]) 
end do: 
phi := subs(y = y-X, phi): 
phi := subs(solve({eval(phi, y = -(1/2)*R+X), subs(y = f, phi)}, {b[0], b[1]}), phi): 
f := piecewise(`and`(z >= -R, z <= 0), z+(1/2)*R+X, -z+(1/2)*R+X): 
Digits := 4: 
2*int(phi*G, [z = y-(1/2)*R .. 0, y = -(1/2)*R .. (1/2)*R], numeric);
Time = time()-T;

@Carl Love 

Thank you very much for your hints

Please consider second branch plot. It seems that boundary conditions are not met.

I change the last rows of the program; as result s(r) at some certain points (N=12) are calculated as follows:

[20., -900],[20.08333333, -838.0971652], [20.16666667, -776.878059], [20.25000000, -716.331788], [20.33333333, -656.447676], [20.41666667, -597.215263], [20.50000000, -538.624300], [20.58333333, -480.664740], [20.66666667, -423.326736], [20.75000000, -366.600637], [20.83333333, -310.476982], [20.91666667, -254.946494], [21.00000000, -200.000079]

I want to check the correctness of the results. How can I validate them by MAPLE?

Boundray conditions : {s(a)=P,s(c)=-200}

@Christian Wolinski 

When I use numeric integration, Maple gives the correct answer.

Perhaps these types of bugs are removed in the new versions.


Thank you for your comments.

The selected region for phi or Integrand was not correct. By rectifying integration domain, it is observed that both the integrand and integral values are posiitive.

int(Integrand*G, [z = y-R .. R-y, y = 0 .. R], numeric);

plots[implicitplot3d](Phi = Integrand, z = -R .. R, y = 0 .. R, Phi = 0 .. 0.1e-1, color = ColorTools[Gradient]("Red" .. "Blue", best)[4], grid = [50, 50, 20]);


Please look over the following code


restart; sub := proc (arg) subs(`$`(d[4-i] = diff(w(r), `$`(r, 4-i)), i = 0 .. 3), d[0] = w(r), arg) end proc; isub := proc (arg) subs(`$`(diff(w(r), `$`(r, 4-i)) = d[4-i], i = 0 .. 3), w(r) = d[0], arg) end proc; df := proc (arg, i) sub(diff(arg, d[i])); (-1)^i*(diff(%, `$`(r, i))) end proc; E := proc (ode) isub(ode); sub(diff(%, d[0]))+df(%, 1)+df(%, 2)+df(%, 3)+df(%, 4) end proc; EL := proc (ode, f) E(ode); f*VectorCalculus:-Laplacian(%/f, 'polar'[r, t]) end proc; eq := proc (ode1, ode2, f) isub(E(ode2)-EL(ode1, f)); factor(%); {seq(coeff(%, d[j]), j = 0 .. 4)} end proc

A := eq(-r*sub(d[1]^2), sub(f[1](r)*d[1]^2+f[2](r)*d[2]^2), r)

{0, 2*(-(diff(f[1](r), r))*r^2-1)/r^2, 2*(2*(diff(f[2](r), r))*r^2-2*r^2)/r^2, 2*(f[2](r)*r^2-r^3)/r^2, 2*((diff(diff(f[2](r), r), r))*r^2-f[1](r)*r^2+r)/r^2}


f[2](r) = solve(A[4], f[2](r))

f[2](r) = r



f[1](r) = 1/r+_C1




Please hint me how it is possible to find an appropriate guess for general form of ode2(r) if ode1(r) is assumed to be q(r)*(w(r))^2, in which q(r) is equal to r*(A+B*r^n)*(1+C*r)^D

D is natural number and A to C and n are real positive constant numerals.



Dear Acer

Thanks alot for your kindness

The method proposed by sand15  is useful.

For example when I use axial rigidity concept (a large amount is assigned to cross sectional area, relatively) the results for some special cases (without drifts) are coincide with the conventional solutions based on the stability functions. This proves the validity of the fsolve results.


Thank you for taking your time.

Your hints are very interesting.

I will use these hints in similar problems (commands like unapply are very useful for large data).


Thanks for your valuable hints. The matrix is symmetric.

Please run the atached code and follow the steps (if you have atleast 10 minutes!):

- Number of the elements = 3

- Number of the nodes = 4

- Elements have the same modulus of elasticity (modulus of elasticity for all elements is 1)

- Elements have the same cross sectional area (cross sectional area for all elements is 1)

- Elements have the same moment of inertia (moment of inertia for all elements is 1)

 X and Y coordinates for nodes are as follows:

node 1 = (0,0)

node 2 = (0,1)

node 3 = (3,4)

node 4 = (4,4)

Node number for begining and end of elements:

element 1 =   begining 1  end 2

element 2 =  begining 2  end 3

element 3 =  begining 3  end 4

number of conditional nodes: 2

conditional node number (for 1 from 2) = 1

select x and y

conditional node number (for 2 from 2) = 4

select x, y, theta

I found the answer with trial and error. In general, how it is possible to get the minimum positive real root without calculating parametric determinant?

For example assume that we use Newton iterative method. The determinant will be calculated at a certain point very fast. But the main problem is that the derivative of the determinant must be calculated at that point symbolically, which is a very time consuming procedure. If it is possible, please propose a way to calculate this part of Newton iterative method (the derivative of determinant at certain value of P in each step) numerically rather than symbolic computing of the determinant derivative.


Thank you for your helpful hints.

@Christian Wolinski

Your idea helps me, but since you reply (you do not use answer option) It is impossible to like your replies (or Thumb them). I do that for your answer, I hope it make you glad.



This code is comprehensive



Thanks for your useful procedure

@Christian Wolinski 




Thank you for taking your time

The answer  4.450852611  is coincide with 15.4182/sqrt(12.)

The small amount of imaginary part may be arisen from somethings like the small number of digits and so on. As you see, the plot command ignore it. I forget to use complex option in fsolve.

Thank you again for your valuable hints


1 2 3 4 5 6 Page 1 of 6