>> From: "madham53"
Hello,
I am solving a system of non-linear equation.
Every time I run the worksheet, I get different answers.
Its an iterative type problem. Though the error here
is only in the second digits, it may blow up when the
iterations are increased.
Is there anything that can be done to get consistent answers
Thank you,
Mahendra
Here's my code:
----------------------------------------------------------------------------
----------
> restart;with(LinearAlgebra):
EVALUATION OF CONSTANTS BEGIN
> m:=1:nu:=0.25:phi:=4:N:=1:alpha:=0:
m = Number of half sine waves in the direction of loading
nu = Poisson's ratio, taken as 0.25
phi = Aspect ratio of plates, taken as 4
N = Compressive load acting at the web-flange intersection producing
only bending
> beta:=(m*Pi/phi)^2:
EVALUATION OF CONSTANTS FOR n = 0
> M[0]:=MatrixInverse(<<0,0,(80-beta*nu),1,(20-beta*nu),(60-5*(2-nu)*beta)
>|<0, 0,
(92-2*beta*nu),2,(12-beta*nu),(24-4*(2-nu)*beta)>|<0,6,(96-4*beta*nu), 4,
(6-beta*nu),(6-3*(2-nu)*beta) >|<
2,0,(64-8*beta*nu),8,(2-beta*nu),-2*(2-nu)*beta>|<0,-(2-nu)*beta,(-16*beta*n
u),16,-beta*nu,-(2-nu)*beta >|<-beta*nu
,0,-32*beta*nu,32,-beta*nu,0>>).Matrix(<<0,0,beta*nu-60,-1/2,beta*nu-30,6*(2
-nu)*beta-120>>):
> evalf(M[0]):
> A[0]:=evalf(M[0][1,1]):
> B[0]:=evalf(M[0][2,1]):
> C[0]:=evalf(M[0][3,1]):
> d[0]:=evalf(M[0][4,1]):
> E[0]:=evalf(M[0][5,1]):
> F[0]:=evalf(M[0][6,1]):
> evalf(M[0]);
[-7.64916722]
[ ]
[14.16520744]
[ ]
[-9.12844676]
[ ]
[1.95534820 ]
[ ]
[-50.7376473]
[ ]
[25.3591287 ]
EVALUATION OF VONKARMAN'S PLATE EQUATIONS
> f[1,0]:= (xi,eta) ->
N/2*(1+alpha-2/3*alpha*eta)*eta^2+((sum(sum('b[p,q]*sin(p*Pi*xi)^2*(eta^(q+6
)-3*eta^(q+5)+13/4*eta^(q+4)-3/2*eta^(q+3)+1/4*eta^(q+2))','q'=0..0),'p'=1..
1))):
> w[1,0]:= (xi,eta) ->
q[0]*sin(Pi*xi)*(eta^6+A[0]*eta^5+B[0]*eta^4+C[0]*eta^3+d[0]*eta^2+E[0]*eta^
1+F[0]*eta^0);
6 5
w[1, 0] := (xi, eta) -> q[0] sin(Pi xi) (eta + A[0] eta
4 3 2
+ B[0] eta + C[0] eta + d[0] eta + E[0] eta + F[0])
CONSTITUENTS OF "A" MATRIX
>
A[1,0][1,0]:=(1/phi^2*diff(f[1,0](xi,eta),[xi,xi,xi,xi])+2*diff(f[1,0](xi,et
a),[xi,xi,eta,eta])+phi^2*diff(f[1,0](xi,eta),[eta,eta,eta,eta]))*(sum(sum('
sin(p*Pi*xi)^2*(eta^(q+6)-3*eta^(q+5)+13/4*eta^(q+4)-3/2*eta^(q+3)+1/4*eta^(
q+2))','q'=0..0),'p'=1..1)):
> Ae[1,0][1,0]:=int(((int(A[1,0][1,0],xi=0..1))),eta=0..1);
Ae[1, 0][1, 0] :=
4 2
1/3843840 b[1, 0] Pi + 1/9240 b[1, 0] Pi + 9/70 b[1, 0]
CONSTITUENTS OF "B" MATRIX
> B[1,0][0,0]:=(-(diff(w[1,0](xi,eta),[xi,eta]))^2 +
diff(w[1,0](xi,eta),[xi,xi])*diff(w[1,0](xi,eta),[eta,eta])) *
(sum(sum('sin(p*Pi*xi)^2*(eta^(q+6)-3*eta^(q+5)+13/4*eta^(q+4)-3/2*eta^(q+3)
+1/4*eta^(q+2))','q'=0..0),'p'=1..1)):
> Be[1,0][0,0]:=int(((int(B[1,0][0,0],xi=0..1))),eta=0..1);
2
Be[1, 0][0, 0] := -3.628759635 q[0]
CALCULATION OF CONSTITUENTS OF "C" MATRIX
CONSTITUENTS OF "C" MATRIX
> C[0,0]:=( 1/phi^2 * diff(w[1,0](xi,eta),[xi,xi,xi,xi])
> +
> 2 * diff(w[1,0](xi,eta),[xi,xi,eta,eta]) +phi^2 * diff(
w[1,0](xi,eta),[eta,eta,eta,eta])) * (sum('(eta^(k+6) + A[k] * eta^(k+5) +
B[k] * eta^(k+4) + C[k] * eta^(k+3) + d[k] * eta^(k+2) + E[k] * eta^(k+1) +
F[k] * eta^k)','k'=0..0)):
> Ce[0,0]:=(int(((int(C[0,0],xi=0..1))),eta=0..1));
Ce[0, 0] := 24826.95903 q[0]
CALCULATION OF CONSTITUENTS OF "D" MATRIX
CONSTITUENTS OF I COLUMN IN "D" MATRIX
> d[1,0][0,0]:= (-12 * (1 - nu^2) * (( diff( f[1,0](xi,eta),[xi,xi] ) *
diff( w[1,0](xi,eta),[eta,eta]) + diff( f[1,0](xi,eta),[eta,eta] ) * diff(
w[1,0](xi,eta),[xi,xi] ) - 2 * diff ( f[1,0](xi,eta),[xi,eta] ) * diff (
w[1,0](xi,eta),[xi,eta] ) ) ) * ( sum ('(eta^(k+6) + A[k] * eta^(k+5) +
B[k] * eta^(k+4) + C[k] * eta^(k+3) + d[k] * eta^(k+2) + E[k] * eta^(k+1) +
F[k] * eta^k)','k'=0..0 ) )):
> de[1,0][0,0]:=(int(((int(d[1,0][0,0],xi=0..1))),eta=0..1));
de[1, 0][0, 0] := 15014.70885 q[0] - 1.159069399 q[0] b[1, 0]
> eqs1:=Ae[1,0][1,0]+Be[1,0][0,0];
4 2
eqs1 := 1/3843840 b[1, 0] Pi + 1/9240 b[1, 0] Pi + 9/70 b[1, 0]
2
- 3.628759635 q[0]
> eqs2:=Ce[0,0]+de[1,0][0,0];
eqs2 := 39841.66788 q[0] - 1.159069399 q[0] b[1, 0]
> solve({eqs1,eqs2});
>
{b[1, 0] = 0., q[0] = 0.},
{q[0] = -35.04662201, b[1, 0] = 34373.84156},
{b[1, 0] = 34373.84156, q[0] = 35.04662201}
----------------------------------------------------------------------------
-----------------
|
(NAG)
