@acer 

restart:

T := time[real]():

 M := 50: 

Digits := 30:

nu:=0.3:   

L := 5*h:  

h := 1:

N := 0.5:

Em := 70:

Ec := 380:  

E :=Em*(1-(y/h+1/2)^N)+Ec*(y/h+1/2)^N:   

R := (1/2)*h: 

G := E/(2*(1+nu)):  

 beta := Pi^2/L^2:   

phi := add(a[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  

a[k+1] := solve(st[k], a[k+1])  

end do:   

bcs := {eval(phi, y = z+1/2), eval(phi, y = 1/2-z)}:   

phi := subs(solve(bcs, {a[0], a[1]}), phi):

inter:=[u = -1/2 .. 1/2, v = -1/2 .. 1/2]:   

intg := 1/2*phi*G:   

intg:=subs(y=(u+v)/2,z=(v-u)/2,intg):   

intg,oldintg:=CodeTools:-Usage(collect(intg,lhs~(inter),uu->simplify(uu,size))),intg:

result := CodeTools:-Usage(evalf(Int(intg,inter,digits=15,epsilon=5e-7,method=_cuhre))):   

result := evalf[7](result);     

Time := time[real]()-T;

@acer 

Dear acer

Thank you very much for your valuable advice

For case III, like other options than defalut, the method you suggested does not work. I change the variables to have numbers as interval for double integration. But it takes too long to get a result. If you have time, please check the code and improve it.

elif cb1=3
then
intg := 1/2*phi*G:
intg:=subs(y=(u+v)/2,z=(v-u)/2,intg):
newintg := collect(intg,[r,t],uu->simplify(uu,size)):
inter:=[u = -1/2 .. 1/2, v = -1/2 .. 1/2]:

@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?

restart;
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.

@tomleslie 

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

Download 01.mw

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.

Thanks

@acer 

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.

@sand15 

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).

@acer 

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.

01.mw

@rlopez 

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.

@mmcdara 

This code is comprehensive

Thanks

@mmcdara 

Thanks for your useful procedure