## 140 Reputation

9 years, 108 days

## The missed value (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;

## 2nd Branch...

Thank you very much for your hints

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

## Check the correctness of calculated amou...

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

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

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

## Thanks...

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]);

## Procedures to facilitate solving of prob...

Please look over the following code

 (1)

 (2)

 (3)

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

## The function definition helps to avoid s...

Dear Acer

The method proposed by 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.

## Thanks...

Thank you for taking your time.

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

## This is Not the eigenvalue problem...

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

## @Christian WolinskiYour idea helps me, b...

@Christian Wolinski

## Thanks again...

This code is comprehensive

Thanks

Thanks

## Thanks...

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
﻿