Unable to solve for a variable

Hi all,

I am trying to solve for a variable 'a' and saving it into an array, However, I am left with an inifinitely long running program. What should I do?

Heres a section of the code:

clc, clear, clear all
J_value=0.5
syms a;
for TonTc=.01:0.01:.99
    sig_sigo_1=((2*J_value+1)/2*J_value)*coth((2*J_value+1)/2*J_value)*a-(1/2*J_value)*coth(a/2*J_value); 
    sig_on_sigo_2=((J_value+1)/3*J_value)*TonTc*a; 
    z=sig_sigo_1-sig_on_sigo_2;
    for k=1:99
        val_a(k)=vpasolve(z,a);
    end
end

Cheers

Answer 1

[Walter Roberson]

J_value=0.5;
syms a;
TonTc_vals = .01:0.01:.99;
num_Tonc = length(TonTc_vals);
val_a = zeros(num_Tonc, 1);

for k = 1 : num_Tonc
    TonTc = TonTc_vals(k);
    sig_sigo_1=((2*J_value+1)/2*J_value)*coth((2*J_value+1)/2*J_value)*a-(1/2*J_value)*coth(a/2*J_value);
    sig_on_sigo_2=((J_value+1)/3*J_value)*TonTc*a;
    z=sig_sigo_1-sig_on_sigo_2;
    sol = vpasolve(z, a, -1);
    if isempty(sol)
        val_a(k) = nan;
    else
        val_a(k) = sol(1);
    end
    if k == 1
        fplot([z, 0], [0.5 1])
        title("plot for TonTc = " + TonTc)
        ylim([-2 2])
    end
end

figure
plot(TonTc_vals, val_a)
title("a vs TonTc")

That first plot makes it clear that there is no root at 0.75

[Jonas Freiheit]

Hi Walter, thanks for the help. For some reason when I try to plot TonTc against sig_sigo I get this plot instead of the theoretical plot that I should obtain. Theoretical attached. Apparantly To plot sig_sigo I need to solve for both 'a' and sig_sigo as separate variables and simply subbing in 'a' doesn't yield the theoretical result for on of the curves of sig_sigo between 0 and 1? What should I do to yield a sig_sigo value between ? Cheers

else

val_a(k) = sol(1);

sig_sigo(k)=((J_value+1)/3*J_value)*TonTc*val_a(k); %This is what I tried but didn't work

[Walter Roberson]

Your denominators were consistently messed up. Everywhere that had something divided by 2*J you had coded as division by 2 and multiplication by J -- x/2*J where you needed x/(2*J)

J_value=0.5;
syms a;
TonTc_vals = .01:0.01:.99;
num_Tonc = length(TonTc_vals);
val_a = zeros(num_Tonc, 1);

for k = 1 : num_Tonc
    TonTc = TonTc_vals(k);
    sig_sigo_1=((2*J_value+1)/(2*J_value))*coth((2*J_value+1)/(2*J_value))*a-(1/(2*J_value))*coth(a/(2*J_value));
    sig_on_sigo_2=((J_value+1)/(3*J_value))*TonTc*a;
    z=sig_sigo_1-sig_on_sigo_2;
    sol = vpasolve(z, a, -1);
    if isempty(sol)
        val_a(k) = nan;
    else
        val_a(k) = sol(1);
    end
    if k == 1
        fplot([z, 0], [0.5 1])
        title("plot for TonTc = " + TonTc)
        ylim([-2 2])
    end
end

figure
plot(TonTc_vals, val_a)
title("a vs TonTc")

[Jonas Freiheit]

Hi Walter,

Really sorry about the late response, I didn't get a notification that there was a response and only checked today. Thanks for the help with that really was a small mistake to not put brackets onto the denominators however, I seem to be getting this problem with my sigma/sigma0 values, I am kinda close but not exact and am unsure what my error in the code is. I just substituted the new 'a' values into the sigma/sigma0 section. Could it be that there are two possible solutions to 'a'?

for k = 1 : num_Tonc
    TonTc = TonTc_vals(k);
    sig_sigo_1=((2*J_value+1)/(2*J_value))*coth((2*J_value+1)/(2*J_value))*a-(1/(2*J_value))*coth(a/(2*J_value));
    sig_on_sigo_2=((J_value+1)/(3*J_value))*(TonTc)*a;
    z=sig_sigo_1-sig_on_sigo_2;
    sol = vpasolve(z, a, -1);
    if isempty(sol)
        val_a(k) = nan;
        sig_sigo(k)=nan;
    else
        val_a(k) = sol(1);
        sig_sigo(k)=(((J_value+1)/(3*J_value))*(TonTc)*val_a(k));

Cheers

[Walter Roberson]

Yes, there are two solutions.

J_value=0.5;
syms a;
TonTc_vals = .01:0.01:.99;
num_Tonc = length(TonTc_vals);
val_a = zeros(num_Tonc, 1);

for k = 1 : num_Tonc
    TonTc = TonTc_vals(k);
    sig_sigo_1=((2*J_value+1)/(2*J_value))*coth((2*J_value+1)/(2*J_value))*a-(1/(2*J_value))*coth(a/(2*J_value));
    sig_on_sigo_2=((J_value+1)/(3*J_value))*TonTc*a;
    z=sig_sigo_1-sig_on_sigo_2;
    sol = vpasolve(z, a, -1);
    if isempty(sol)
        val_a(k) = nan;
    else
        val_a(k) = sol(1);
    end
    if k == 1
        fplot([z, 0], [-2 2])
        title("plot for TonTc = " + TonTc)
        ylim([-2 2])
    end
end

figure
plot(TonTc_vals, val_a)
title("a vs TonTc")

[Jonas Freiheit]

Yes, Sorry How would I plug in the second solutions?

Just trying to test to see if that gives me the theoretical results. This is the equation I'm supposed to follow to obtain sigma/sigma0 which I think I've correctly used in the code.

Answer 2

[David Hill]

What are you solving? You need to set z equal to something.

z=sig_sigo_1-sig_on_sigo_2==0;%what do you want z to be?

[John D'Errico]

Actually, solve implicitly assumes zero if you specify nothing. For example

syms X
solve(X-1,X)

ans = 

1

So solve implicitly assume the equation was x-1 == 0, despite the lack of any right hand side provided.

[David Hill]

Thanks, good to know.

[Jonas Freiheit]

Hi,

For some reason when I do z==0 and do vpasolve(z,a) for TonTc=0.1 I get 0.972 when the actual correct answer is 0.75 from my calculator? Howcome is this wrong?

Thanks

Also solve doesn't work and has to be vpasolve