【求助】matlab求解非线性方程组问题。
本帖最后由 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=;
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=;
=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 ()
ylabel('Saturation (-)','Fontsize',12);
xlabel('RH (-)','Fontsize',12);
box on
关于这个问题，可能的原因有两个，
1. 程序书写错误。
2. 方程组无解，也就是说不存在一条抛物线既经过A点又与黑线相切。
如果是第二个原因，那么用什么方法可以证明抛物线不可能和黑线相切（在内侧相切）？
页:
[1]