Labfans是一个针对大学生、工程师和科研工作者的技术社区。 论坛首页 | 联系我们(Contact Us)
MATLAB爱好者论坛-LabFans.com
返回   MATLAB爱好者论坛-LabFans.com > 工程数学软件 > MATLAB论坛
MATLAB论坛 一切MATLAB相关问题在此讨论。
回复
 
主题工具 显示模式
旧 2007-09-01, 19:40   #1
小妮妮
初级会员
 
注册日期: 2007-09-01
帖子: 1
声望力: 0
小妮妮 正向着好的方向发展
默认 【求助】超越方程组的fsolve解法

我用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);
上传的附件
文件类型: doc 公式.doc (52.0 KB, 25 次查看)
小妮妮 当前离线   回复时引用此帖
旧 2008-07-11, 15:33   #2
kissiopo
初级会员
 
注册日期: 2008-01-09
年龄: 41
帖子: 2
声望力: 0
kissiopo 正向着好的方向发展
默认

我现在也蛮烦躁 在想么样求解比较好 非线性方程组 那个误差对于初值不同 偏差蛮大 :cry:
kissiopo 当前离线   回复时引用此帖
旧 2008-07-16, 15:03   #3
mmfang
初级会员
 
注册日期: 2008-07-16
年龄: 41
帖子: 1
声望力: 0
mmfang 正向着好的方向发展
默认

帮顶 我也是被非线性方程组卡着了
mmfang 当前离线   回复时引用此帖
旧 2008-09-21, 20:11   #4
未注册
游客
 
帖子: n/a
默认 回复: 【求助】超越方程组的fsolve解法

同志啊!!!!!!
  回复时引用此帖
回复

主题工具
显示模式

发帖规则
不可以发表新主题
不可以发表回复
不可以上传附件
不可以编辑自己的帖子

启用 BB 代码
论坛启用 表情符号
论坛启用 [IMG] 代码
论坛禁用 HTML 代码



所有时间均为北京时间。现在的时间是 17:36


Powered by vBulletin
版权所有 ©2000 - 2025,Jelsoft Enterprises Ltd.