Answer not correct when comparing to wolframalpha

Hello,

I've never had this problem before, but the last thing preventing me from using a program I wrote is that the term 1 that I've isolated gives the wrong answer for some unknown reason. It's line 56 where I just made the variable "term 1" used in line 54 in the real formula. I put in the correct number at the end of the Function_M into wolframalpha to find the correct plot and it just seems to not give the correct answer when function = 0.

clear;
clc;
%Input Values
%Shear_max=560*(10^6); %yield stress from data sheet
g_max_A=172.369*(10^6); %max cyclic shear stress austentite MPa (25ksi)
g_max_M=68.9476*(10^6); %max cyclic shear stress martensite MPa (10ksi)
G_M=10.8*(10^9); %Shear Modulus of the cold Material
G_A=28.8*(10^9); %Shear Modulus of the hot Material
v=0.33; %Poisson's ratio
%Wire and Coil Diameter
d_max=0.058444; %Stroke Displacement (increased by 1/(%usable strain))
d=150*(10^-6); %The diameter of the Wire (microns converted to m) cools in 1 second with mesh geometry
F_M=0.9807/2; %100 grams in newtons / 2 for 50 grams martensite pulling force
F_A=0.9807/2*(g_max_A/g_max_M); %Force adjusted for austenite max stress at same D
%w=[1 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30];
w=10;
D=g_max_M*pi*d^3.*w./(8*F_M); %Spring diameter using martensite numbers
%D_A=g_max_A*pi*d^3.*w./(8*F_A);
D_Mandrel=D-d;
g_eight=0.71;
g_twelve=1;
d_net=0;
c=1;
length_f=0.12;
        angle_A_f=sym('angle_A_f');
        angle_M_f=sym('angle_M_f');
        Function_A=sym('Function_A');
        Function_M=sym('Function_M');
%Iterating final length to find a diameter for force
%while (length_f >= 0.120)
%    length_f=length_f-0.001;
angle_i=40;
while length_f >= 0.100
        angle_i=angle_i-1; %iterated initial angle
        angle_A_f=sym('angle_A_f');
        angle_M_f=sym('angle_M_f');
%        Function_A=0;
        Function_A=G_A*d^4*pi/(8*D^2)*...
             cosd(angle_i)^2*...
             (sind(angle_A_f)-sind(angle_i))/...
             (cosd(angle_A_f)^2*...
             (cosd(angle_A_f)^2+sind(angle_A_f)^2/...
             (1+v)))-...
             (F_A/w);
        %angle_A_f=double(angle_A_f);
%        if angle_i==1
        angle_A_f=double(vpasolve(Function_A==0, angle_A_f, [0 90]));
%        else
        %angle_A_f=double(vpasolve(Function_A==0, angle_A_f, [0 90]));
%        end
        Function_M=G_M*d^4*pi/(8*D^2).*cosd(angle_i)^2*(sind(angle_M_f)-sind(angle_i))/...
             (cosd(angle_M_f)^2*...
             (cosd(angle_M_f)^2+(sind(angle_M_f)^2/...
             (1+v))))-(pi*d^3/(8*D)*G_M*0.06*g_eight) -...
             (F_M/w);
        term1=-(G_M.*g_eight);
        term1_2=(pi*d^3)/(8*D);
        term2=-(F_M/w);
        angle_M_f=double(vpasolve(Function_M==0, angle_M_f, [0 90]));
        
        n=cosd(angle_i)*d_max/(pi*D.*(sind(angle_M_f)-sind(angle_A_f)));
        length_f=n*pi*D*(tand(angle_M_f));
        length_i=n*pi*D*(tand(angle_i));
        %pitch_f=tand(angle_f)*pi*D; %Coil pitch
         %Length of unsprung spring set to forearm
        %length_i=length_f-dl; %Final length = initial + displacement
        %n=length_f/pitch_f; %Number of loops in coil
        %angle_i=atand(length_i/(pi*n*D));
        %Shear Strain and Detwinned Martensite Fraction
%        angle_A_f0 = 0;
%        for c=1:length(w)
%        angle_A_f = fzero(@(angle_A_f) findangle_A(angle_A_f,G_A,d,D_A,angle_i,v,w,c)...
%            ,angle_A_f0);
%        
%        angle_M_f0 = 0;
%        angle_M_f = fzero(@(angle_M_f) findangle_M(angle_M_f,G_M,d,D_M,angle_i,v,g_eight,w,c)...
%            ,angle_M_f0);
%        for c=1:length(w)

%        end
        if angle_i <= 3
            break
        end
end

Here are some pictures of the formula used on line 51 and proof from wolframalpha that the numbers are different.

Formula at angle_i=5:

0=10.8^9*(0.15/1000)^4/(8*0.00186^2)*(cos(5 degrees)^2)*(sin(x degrees)-(sin(5 degrees)))/(sin(x degrees)^2*(cos(x degrees)^2+sin(x degrees)^2/(1+0.33)))-0.3762

Formula that I'm trying to reproduce:

everything besides yL and ESt is referenced clearly in the code, and yL = 0.06, and ESt = 0.71.

[Alan Stevens]

Your WolframAlpha input doesn't match your code, nor the mathematical formula. For example,

10.8^9 (WolframAlpha) is not the same as 10.8*(10^9) (your code)

You seem to have dropped a pi multiplying the first term in WolframAlpha.

[Vincent Pipitone]

Ok, I fixed the 10^9 term and added the pi to wolframalpha. The term2 at the end was fixed to give the correct value of -0.3762. The graphs still do not match to suggest the answer is correct. I am getting 30.7266 degrees on matlab and 5.299 degrees on wolframalpha. Here are the images of the formulas and the formula text:

0=10.8*10^9*pi*(0.15/1000)^4/(8*0.00186^2)*(cos(5 degrees)^2)*(sin(x degrees)-(sin(5 degrees)))/(sin(x degrees)^2*(cos(x degrees)^2+sin(x degrees)^2/(1+0.33)))-0.3762

Thank you for your help so far!

[Alan Stevens]

Well, if I copy what you've put in WolframAlpha to MATLAB, like so

f = @(x) 10.8*10^9*pi*(0.15/1000)^4/(8*0.00186^2)*(cosd(5)^2)*(sind(x) ...
 -(sind(5)))./(sind(x).^2.*(cosd(x).^2+sind(x).^2/(1+0.33)))-0.3762;

x = -210:210;
y = f(x);
plot(x,y), grid
axis([-210 210 -4 1.5])

I get the following

This looks the same as the WolframAlpha graph to me! There must be some other difference in your original Matlab code - I'll have a closer look.

Answer 1

[Alan Stevens]

Your WolframAlpha expression is a mixture of your Function_A and Function_M. You have used G_M, from Function_M, but in the denominator you have

sin(x degrees)^2*(cos(x degrees)^2+sin(x degrees)^2/(1+0.33))

which is from Function_A.

Function_M has

cos(x degrees)^2*(cos(x degrees)^2+sin(x degrees)^2/(1+0.33))

[Vincent Pipitone]

Thank you very much! I am getting the same answer in wolframalpha as I am in the program, which suggests everything is typed in correctly.

[John D'Errico]

Alan has done a great deal of work here, tracking down the problem. I hope you will accept his answer, and provide an upvote.