马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
本帖最后由 gulley 于 2012-4-26 20:36 编辑
大家好。初学matlab,试图用它求解一组非线性方程组,总是得不到合理的结果。坛子里有人推荐1stopt,我也试了,每次的结果都不一样。
这里,我先说明一下我要做什么,再附上程序,希望大家指导一下。
1. 求解问题
我有两条曲线(图中蓝线和黑线),我准备从一条曲线上(蓝线)的一点(点A)做一条抛物线,使这条抛物线与另外一条曲线(黑线)相切。而在起点处,我控制这条抛物线的斜率。
所以,共有三条曲线,他们的公式为(除x,y外,其他都是常数):
从上面的说明可以看出,已知起点(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
|