声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 1599|回复: 1

[综合讨论] 【求助】matlab求解非线性方程组问题。

[复制链接]
发表于 2012-4-26 17:50 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
本帖最后由 gulley 于 2012-4-26 20:36 编辑

大家好。初学matlab,试图用它求解一组非线性方程组,总是得不到合理的结果。坛子里有人推荐1stopt,我也试了,每次的结果都不一样。
这里,我先说明一下我要做什么,再附上程序,希望大家指导一下。

1. 求解问题
我有两条曲线(图中蓝线和黑线),我准备从一条曲线上(蓝线)的一点(点A)做一条抛物线,使这条抛物线与另外一条曲线(黑线)相切。而在起点处,我控制这条抛物线的斜率。
所以,共有三条曲线,他们的公式为(除x,y外,其他都是常数):


捕获.JPG

从上面的说明可以看出,已知起点(x1,y1)和起点处斜率(Ks)。未知 A,B,C和终点(只需知道x3即可)。我们有四个条件,可以列四个方程,解出这个四个未知数。分别为:
(1)起点在抛物线上:y1=A*x1.^2.+B*x1+C
(2)起点处给定斜率:2*A*x1+B = Ks
(3)终点在抛物线上,也在黑线上:A*x3.^2.+B*x3+C = ((para1/alpha_s*log(x3)).^(beta_s/(beta_s-1.))+1.).^(-1./beta_s)
(4)终点处斜率与黑线的斜率相同:2*A*x3+B = diff(((para1/alpha_s*log(x3)).^(beta_s/(beta_s-1.))+1.).^(-1./beta_s))

2. 程序代码
以下程序中 rh=x,s=y。
起点为 (0.3,0.3),起点斜率为 Ks=0.2。
% 以下代码保存在scanning.m文件中
function f = scanning(x,Ks,rh1,s1)

syms A B C rh3 rh rh_cross kk;
global kk para1 alpha_s beta_s alpha_e beta_e;
A=x(1);
B=x(2);
C=x(3);
rh3=x(4);

y1=diff(rh);
y2=diff(((para1/alpha_e*log(rh)).^(beta_e/(beta_e-1.))+1.).^(-1./beta_e)); % 求带有变量的函数的导数

f1=A*rh1.^2.+B*rh1+C-s1; % 起点在抛物线上
kk=Ks*subs (y1,rh,rh1);
f2=2*A*rh1+B-Ks*subs(y1,rh,rh1); %起点处的斜率
f2=2*A*rh1+B-kk; %抛物线起点处的斜率
f3=A*rh3*rh3+B*rh3+C-(((para1/alpha_s*log(rh3)).^(beta_s/(beta_s-1.))+1.).^(-1./beta_s)); %终点在抛物线上,并且在主曲线上
f4=2*A*rh3+B-subs(y2,rh,rh3); %抛物线在终点处的斜率与主曲线相同

f=[f1;f2;f3;f4];

end

%%%%%%%%%%%%%%%%%%
% 以下代码另存一个文件
clear all
global kk para1 alpha_s beta_s alpha_e beta_e;

Ks = 0.2;   % initial slope
para1= -1.3694e+008;
alpha_s =32.70e6;
beta_s = 2.16; %
alpha_e = 1.48e6;
beta_e = 3.79;
rh1 = 0.3;  % 给定起始点
s1 = rh1

%solve the nonlinear eqautions using an extra function
options=optimset('Display','iter');
x0=[1.0;1.0;1.0;0.5];
[x,fval]=fsolve(@(x) Bjorn_scanning(x,Ks,rh1,s1),x0,options)

% give the result for each parameter
A=x(1);
B=x(2);
C=x(3);
rh3=x(4);

%calculate the curve
rh_sc=(rh1:0.001:rh3);
s_sc=A*rh_sc.^2+B*rh_sc+C;

%calculate main curves
rh_st=(0.01:0.01:1.);
s_st=rh_st;
s_en=((para1/alpha_e*log(rh_st)).^(beta_e/(beta_e-1.))+1.).^(-1./beta_e);

%
figure(1)
hold on;
plot(rh_sc,s_sc,'r');
plot(rh_st,s_st,'b');
plot(rh_st,s_en,'k');
axis ([0 1 0 1])
ylabel('Saturation (-)','Fontsize',12);
xlabel('RH (-)','Fontsize',12);
box on
untitled.jpg





回复
分享到:

使用道具 举报

 楼主| 发表于 2012-4-26 20:34 | 显示全部楼层
关于这个问题,可能的原因有两个,
1. 程序书写错误。
2. 方程组无解,也就是说不存在一条抛物线既经过A点又与黑线相切。

如果是第二个原因,那么用什么方法可以证明抛物线不可能和黑线相切(在内侧相切)?
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-13 14:00 , Processed in 0.139587 second(s), 22 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表