Unable to solve for a variable

2 次查看(过去 30 天)
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

采纳的回答

Walter Roberson
Walter Roberson 2021-10-7
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
  5 个评论
Walter Roberson
Walter Roberson 2021-10-10
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
Jonas Freiheit 2021-10-10
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.

请先登录,再进行评论。

更多回答(1 个)

David Hill
David Hill 2021-10-6
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?
  3 个评论
Jonas Freiheit
Jonas Freiheit 2021-10-7
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

请先登录,再进行评论。

产品

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by