E1=181e9;E2=10.3e9;G12=7.17e9;s11=12e6;s22=-46e6;t12=-19e6;sxy=[s11;s22;t12];v12=0.28;v21=(E2/E1)*v12;S11=1/E1;S12=-v12/E1;S21=-v21/E2;S22=1/E2;S33=1/G12;S0=[S11 S12 0;S21 S22 0;0 0 S33];C=inv(S0);% T is transform matrix%
for theta=[-90:15:90]m=cosd(theta);n=sind(theta);R=[1 0 0;0 1 0;0 0 2];T=[m^2 n^2 2*m*n;n^2 m^2 -2*m*n;-m*n n*m (m^2-n^2)] ;ST=R*T*S0*(T^-1)*(R^-1)GT=inv(ST)Qxy=inv(T)*inv(S0)*R*T*inv(R)% Ex, Ey, nu_xy, Gxy,Axy %
Sxy=inv(Qxy);Ex=1/Sxy(1,1)Ey=1/Sxy(2,2)Gxy=1/Sxy(3,3)vxy=-Sxy(1,2)*Exend
MATLAB: Store the results(Ex, Ey, vxy) for each iteration in the loop. So that I can plot later
arrayMATLABplottingstore
Related Solutions
Not all nasty looking, high order nonlinear systems of equations have solutions. vpasolve is a capable numerical solver, but it has limits. One issue is that when working in high precision, things take longer to do. All computations are slower. It will probably happen that vpasolve will return a high precision solution eventually.
So I tried posing this as a problem for fsolve. Working in double precision is far faster. fsolve can sometimes have issues with complex valued functions. But it often does seem to work.
Here is your code, pasted into a function to be used for fsolve.
function eqs = funs(S)S11 = S(1);S22 = S(2);S33 = S(3);S12 = S(4);S13 = S(5);S23 = S(6);A11= 0.2346 + 0.1345i;A21= 0.8550 - 0.4363i;A22= -0.0412 + 0.2931i;B11= 0.2346 + 0.1345i;B21= 0.8550 - 0.4363i;B22= -0.0412 + 0.2931i;D11= 0.2346 + 0.1345i;D21= 0.8550 - 0.4363i;D22= -0.0412 + 0.2931i;c11= 0.169033453278504 + 0.127942837737115i;c22= 0.169951074114244 + 0.148423165394202i;c33= 0.168939745364374 + 0.142538799523638i;c21= -0.00172796568379862 + 0.201091742020021i;c31= -0.00152297772077853 + 0.199585730377658i;c32= -0.00523160858486265 + 0.199472640309087i;eqs = [(A11.*A22.*S11 - A21.^2.*S11 - A11 + A11.*B22.*S22 + A11.*D22.*S33 - A21.^2.*B22.*S12.^2 - A21.^2.*D22.*S13.^2 + A11.*A22.*B22.*S12.^2 + A11.*A22.*D22.*S13.^2 + A11.*B22.*D22.*S23.^2 + A21.^2.*B22.*S11.*S22 + A21.^2.*D22.*S11.*S33 + A21.^2.*B22.*D22.*S11.*S23.^2 + A21.^2.*B22.*D22.*S13.^2.*S22 + A21.^2.*B22.*D22.*S12.^2.*S33 - A11.*A22.*B22.*S11.*S22 - A11.*A22.*D22.*S11.*S33 - A11.*B22.*D22.*S22.*S33 - A11.*A22.*B22.*D22.*S11.*S23.^2 - A11.*A22.*B22.*D22.*S13.^2.*S22 - A11.*A22.*B22.*D22.*S12.^2.*S33 - 2.*A21.^2.*B22.*D22.*S12.*S13.*S23 - A21.^2.*B22.*D22.*S11.*S22.*S33 + 2.*A11.*A22.*B22.*D22.*S12.*S13.*S23 + A11.*A22.*B22.*D22.*S11.*S22.*S33)./(A22.*S11 + B22.*S22 + D22.*S33 + A22.*B22.*S12.^2 + A22.*D22.*S13.^2 + B22.*D22.*S23.^2 - A22.*B22.*S11.*S22 - A22.*D22.*S11.*S33 - B22.*D22.*S22.*S33 - A22.*B22.*D22.*S11.*S23.^2 - A22.*B22.*D22.*S13.^2.*S22 - A22.*B22.*D22.*S12.^2.*S33 + 2.*A22.*B22.*D22.*S12.*S13.*S23 + A22.*B22.*D22.*S11.*S22.*S33 - 1)-c11; -(A21.*B21.*(S12 + D22.*S13.*S23 - D22.*S12.*S33))./(A22.*S11 + B22.*S22 + D22.*S33 + A22.*B22.*S12.^2 + A22.*D22.*S13.^2 + B22.*D22.*S23.^2 - A22.*B22.*S11.*S22 - A22.*D22.*S11.*S33 - B22.*D22.*S22.*S33 - A22.*B22.*D22.*S11.*S23.^2 - A22.*B22.*D22.*S13.^2.*S22 - A22.*B22.*D22.*S12.^2.*S33 + 2.*A22.*B22.*D22.*S12.*S13.*S23 + A22.*B22.*D22.*S11.*S22.*S33 - 1)-c21; -(A21.*D21.*(S13 + B22.*S12.*S23 - B22.*S13.*S22))./(A22.*S11 + B22.*S22 + D22.*S33 + A22.*B22.*S12.^2 + A22.*D22.*S13.^2 + B22.*D22.*S23.^2 - A22.*B22.*S11.*S22 - A22.*D22.*S11.*S33 - B22.*D22.*S22.*S33 - A22.*B22.*D22.*S11.*S23.^2 - A22.*B22.*D22.*S13.^2.*S22 - A22.*B22.*D22.*S12.^2.*S33 + 2.*A22.*B22.*D22.*S12.*S13.*S23 + A22.*B22.*D22.*S11.*S22.*S33 - 1)-c31; (A22.*B11.*S11 - B21.^2.*S22 - B11 + B11.*B22.*S22 + B11.*D22.*S33 - A22.*B21.^2.*S12.^2 - B21.^2.*D22.*S23.^2 + A22.*B11.*B22.*S12.^2 + A22.*B11.*D22.*S13.^2 + B11.*B22.*D22.*S23.^2 + A22.*B21.^2.*S11.*S22 + B21.^2.*D22.*S22.*S33 + A22.*B21.^2.*D22.*S11.*S23.^2 + A22.*B21.^2.*D22.*S13.^2.*S22 + A22.*B21.^2.*D22.*S12.^2.*S33 - A22.*B11.*B22.*S11.*S22 - A22.*B11.*D22.*S11.*S33 - B11.*B22.*D22.*S22.*S33 - A22.*B11.*B22.*D22.*S11.*S23.^2 - A22.*B11.*B22.*D22.*S13.^2.*S22 - A22.*B11.*B22.*D22.*S12.^2.*S33 - 2.*A22.*B21.^2.*D22.*S12.*S13.*S23 - A22.*B21.^2.*D22.*S11.*S22.*S33 + 2.*A22.*B11.*B22.*D22.*S12.*S13.*S23 + A22.*B11.*B22.*D22.*S11.*S22.*S33)./(A22.*S11 + B22.*S22 + D22.*S33 + A22.*B22.*S12.^2 + A22.*D22.*S13.^2 + B22.*D22.*S23.^2 - A22.*B22.*S11.*S22 - A22.*D22.*S11.*S33 - B22.*D22.*S22.*S33 - A22.*B22.*D22.*S11.*S23.^2 - A22.*B22.*D22.*S13.^2.*S22 - A22.*B22.*D22.*S12.^2.*S33 + 2.*A22.*B22.*D22.*S12.*S13.*S23 + A22.*B22.*D22.*S11.*S22.*S33 - 1)-c22; -(B21.*D21.*(S23 + A22.*S12.*S13 - A22.*S11.*S23))./(A22.*S11 + B22.*S22 + D22.*S33 + A22.*B22.*S12.^2 + A22.*D22.*S13.^2 + B22.*D22.*S23.^2 - A22.*B22.*S11.*S22 - A22.*D22.*S11.*S33 - B22.*D22.*S22.*S33 - A22.*B22.*D22.*S11.*S23.^2 - A22.*B22.*D22.*S13.^2.*S22 - A22.*B22.*D22.*S12.^2.*S33 + 2.*A22.*B22.*D22.*S12.*S13.*S23 + A22.*B22.*D22.*S11.*S22.*S33 - 1)-c32; (A22.*D11.*S11 - D21.^2.*S33 - D11 + B22.*D11.*S22 + D11.*D22.*S33 - A22.*D21.^2.*S13.^2 - B22.*D21.^2.*S23.^2 + A22.*B22.*D11.*S12.^2 + A22.*D11.*D22.*S13.^2 + B22.*D11.*D22.*S23.^2 + A22.*D21.^2.*S11.*S33 + B22.*D21.^2.*S22.*S33 + A22.*B22.*D21.^2.*S11.*S23.^2 + A22.*B22.*D21.^2.*S13.^2.*S22 + A22.*B22.*D21.^2.*S12.^2.*S33 - A22.*B22.*D11.*S11.*S22 - A22.*D11.*D22.*S11.*S33 - B22.*D11.*D22.*S22.*S33 - A22.*B22.*D11.*D22.*S11.*S23.^2 - A22.*B22.*D11.*D22.*S13.^2.*S22 - A22.*B22.*D11.*D22.*S12.^2.*S33 - 2.*A22.*B22.*D21.^2.*S12.*S13.*S23 - A22.*B22.*D21.^2.*S11.*S22.*S33 + 2.*A22.*B22.*D11.*D22.*S12.*S13.*S23 + A22.*B22.*D11.*D22.*S11.*S22.*S33)./(A22.*S11 + B22.*S22 + D22.*S33 + A22.*B22.*S12.^2 + A22.*D22.*S13.^2 + B22.*D22.*S23.^2 - A22.*B22.*S11.*S22 - A22.*D22.*S11.*S33 - B22.*D22.*S22.*S33 - A22.*B22.*D22.*S11.*S23.^2 - A22.*B22.*D22.*S13.^2.*S22 - A22.*B22.*D22.*S12.^2.*S33 + 2.*A22.*B22.*D22.*S12.*S13.*S23 + A22.*B22.*D22.*S11.*S22.*S33 - 1)-c33];
After saving this function on my search path, I used fsolve, as follows:
format long gSfinal = fsolve(@funs,(1+i)*ones(1,6))Equation solved.fsolve completed because the vector of function values is near zeroas measured by the default value of the function tolerance, andthe problem appears regular as measured by the gradient.<stopping criteria details>Sfinal =Columns 1 through 3 -0.0576463452402932 - 0.0753664180207837i -0.0744682052419572 - 0.0630375797716533i -0.0698118920703255 - 0.0672129482711116iColumns 4 through 6 -0.185173231731065 + 0.109260636702373i -0.183468124937688 + 0.10862578397362i -0.186854054845906 + 0.104967752613526i
As a test that this is a solution in double precision, evaluate funs at that point.
funs(Sfinal)ans = 2.04330996567137e-12 + 8.08547673258886e-13i 2.06829562157673e-12 + 8.28642710004601e-13i 2.07329856408145e-12 + 8.41465785939022e-13i 2.07303618715571e-12 + 7.6394446324457e-13i 2.08355208086708e-12 + 8.14071032806396e-13i 2.05693795329864e-12 + 7.88147325181399e-13i
Note that when calling fsolve, it was necessary to use a complex start point. Even then, I have seen cases where fsolve fails for complex valued problems, though I have not yet pinned down why it may fail. It did seem to work here though.
theta=[-90:15:90];N=length(theta);[...]for i=1:N m=cosd(theta);n=sind(theta);
why bother iterating over the length of theta when you are going to operate on all of theta inside the loop?
Best Answer