小妮妮
2007-09-01, 19:40
我用matlab的fsolve函数求解三元超越方程组，但是得到的解和文献上不同，困扰几个月了，在其他的论坛也求助了，但是一直都找不到解决的办法，现在课题全卡住了，现将代码和公式附上，请高手帮忙看看,万分感谢！sx1,sx2,tc是未知数，在程序运行时部分解输出的exitflag<0。
%%%%%%%%%%主程序%%%%%%%%%%%%%%%%%%%%%%
clear all;clc;close all;
global js;
global q;
Tc =[];
Sx1=[];
Sx2=[];
%q_par = [0.9];
%q_par = [1.065];
q_par = [2.0];
js_par = 2.0:-0.1:0;
n0= [ unifrnd(0,0.4) unifrnd(0,0.4) 1.5];
for q = q_par
js = js_par(1);
options=optimset('display','on');
[n(1,:),fval,exitflag]=fsolve(@myfun,n0,options)
for i = 2:length(js_par)
js = js_par(i);
[n(i,:),fval,exitflag]=fsolve(@myfun,n(i-1,:),options)
end
Tc = [Tc n(:,3)]
Sx1=[Sx1 n(:,1)]
Sx2=[Sx2 n(:,2)]
end
plot(Tc,js_par)
str =[];
for i=1:length(q_par)
str = strvcat(str,num2str(q_par(i)));
end
legend(str,4);
xlabel('Tc/Ja')
ylabel('Js/Ja')
%%%%%%%%%%%%%函数%%%%%%%%%%%%%%%%%%%
function y=myfun(n)
global js
global q
syms tc sx1 sx2
je=1.0;
qs=2.0;
k=1.0;
sx1=n(1);
sx2=n(2);
tc=n(3);
w31=sqrt(qs^2+4*js*qs*sx1+je*q*sx2);
w32=sqrt(q^2+4*q*sx2+je*qs*sx1+je*qs*sx1);
y(1)=qs/(2*w31)*tanh(w31/(2*k*tc))-sx1;
y(2)=q/(2*w32)*tanh(w32/(2*k*tc))-sx2;
N=[2*w31*coth(w31/(2*k*tc))-4*js, -je, 0;
-je, 2*w32*coth(w32/(2*k*tc))-4, -je;
0, -je, 2*w31*coth(w31/(2*k*tc))-4*js];
y(3)=det(N);
%%%%%%%%%%主程序%%%%%%%%%%%%%%%%%%%%%%
clear all;clc;close all;
global js;
global q;
Tc =[];
Sx1=[];
Sx2=[];
%q_par = [0.9];
%q_par = [1.065];
q_par = [2.0];
js_par = 2.0:-0.1:0;
n0= [ unifrnd(0,0.4) unifrnd(0,0.4) 1.5];
for q = q_par
js = js_par(1);
options=optimset('display','on');
[n(1,:),fval,exitflag]=fsolve(@myfun,n0,options)
for i = 2:length(js_par)
js = js_par(i);
[n(i,:),fval,exitflag]=fsolve(@myfun,n(i-1,:),options)
end
Tc = [Tc n(:,3)]
Sx1=[Sx1 n(:,1)]
Sx2=[Sx2 n(:,2)]
end
plot(Tc,js_par)
str =[];
for i=1:length(q_par)
str = strvcat(str,num2str(q_par(i)));
end
legend(str,4);
xlabel('Tc/Ja')
ylabel('Js/Ja')
%%%%%%%%%%%%%函数%%%%%%%%%%%%%%%%%%%
function y=myfun(n)
global js
global q
syms tc sx1 sx2
je=1.0;
qs=2.0;
k=1.0;
sx1=n(1);
sx2=n(2);
tc=n(3);
w31=sqrt(qs^2+4*js*qs*sx1+je*q*sx2);
w32=sqrt(q^2+4*q*sx2+je*qs*sx1+je*qs*sx1);
y(1)=qs/(2*w31)*tanh(w31/(2*k*tc))-sx1;
y(2)=q/(2*w32)*tanh(w32/(2*k*tc))-sx2;
N=[2*w31*coth(w31/(2*k*tc))-4*js, -je, 0;
-je, 2*w32*coth(w32/(2*k*tc))-4, -je;
0, -je, 2*w31*coth(w31/(2*k*tc))-4*js];
y(3)=det(N);