帮忙matlab求两条曲线交点程序,不知问题出在哪里.菜鸟+新手弄不明白.% temperature diffusion 温度扩散曲线z = 0:300; %depth in mbsf 深度t1 = 500; % timeb = 0.054; %thermal gradientDT = 2; %temperature changek = 3E-7; %the

帮忙matlab求两条曲线交点程序,不知问题出在哪里.菜鸟+新手弄不明白.% temperature diffusion 温度扩散曲线z = 0:300; %depth in mbsf 深度t1 = 500; % timeb = 0.054; %thermal gradientDT = 2; %temperature changek = 3E-7; %the
帮忙matlab求两条曲线交点程序,不知问题出在哪里.
菜鸟+新手弄不明白.
% temperature diffusion 温度扩散曲线
z = 0:300; %depth in mbsf 深度
t1 = 500; % time
b = 0.054; %thermal gradient
DT = 2; %temperature change
k = 3E-7; %thermal diffusivity
To = 2.2; %initial temperature
T1 = To+b.*z+DT.*erfc(z./2.*sqrt(k.*t1)); % 温度随深度扩散曲线
hold on;
%plot temperature diffusion line
plot(T1,z,'g'); %green
%Xlabel('Temperature (C)');
%Ylabel('Depth (mbsf)');
% hydrate stability Line 稳定方程曲线
% 3.35% salt water 100% pure methane
P = 0.101325+1030.*9.81.*(750+z).*0.000001; %净水压力
LP = log10(P);
TB = 1./(0.00379-0.000283.*LP)-273.15; %稳定方程曲线
% plot hydrate equilibrium
plot(TB,z,'k--'); %black dash
%到这里都没问题
% junction point T z 求两条曲线交点,下面这行出错.
result=solve('T = To+b.*z+DT.*erfc(z./2.*sqrt(k.*t1))','T = 1./(0.00379-0.000283.*log10(0.101325+1030.*9.81.*(750+z).*0.000001))-273.15');
x = double(result,T);
Y = double(result,z);

帮忙matlab求两条曲线交点程序,不知问题出在哪里.菜鸟+新手弄不明白.% temperature diffusion 温度扩散曲线z = 0:300; %depth in mbsf 深度t1 = 500; % timeb = 0.054; %thermal gradientDT = 2; %temperature changek = 3E-7; %the

您好,

1. 如果您想用solve算子来得出交点值,您需要把您方程适中除自变量外所有的其它参数用数字的形式直接写出来,而不是继续用变量名!

请把您在注释% junction point T z 之后的语句替换为下面的语句 即可成功运行：

zR = solve('1/(0.00379-0.000283*log10(0.101325+1030*9.81*(750+z)*0.000001))-273.15-(2.2+0.054*z+2*erfc(z/2*sqrt(3E-7*500)))','z');

tR = 1/(0.00379-0.000283*log10(0.101325+1030*9.81*(750+zR)*0.000001))-273.15;

zSol = double(zR);

tSol = double(tR);

str = sprintf('The intersection of 2 lines is : \n z = %f \n t = %f', zSol,tSol);

disp(str);

plot(tSol,zSol,'r+','LineWidth',3,'MarkerSize',13);

P.S. 如果您需要查看十分十分确切的结果.请直接在命令窗口(Command Window)中输入：zR 或者tR,便会显示其小数点后二十多位的数值结果.

        zR  = 162.4976393478409132553141680058

        tR  = 11.293566711605738946287734914997

2.我看了您上面写的内容.感觉您应该是属于工程领域的吧?由于,我是搞工科的.所以这样的话,我觉得您并不需要使用solve算子来求出理论交点.即使求出来了,也没有太大的实际意义.往往我们需要的数值顶多在小数点后3~6位左右就差不多了.太深的数值精度对于工科而言毫无意义.

3.基于上面的第二点,我修改了您的命令文档.然后通过赋予容忍度1e-5来确定交点.主要思想是：先重新定义z向量（在这里我是用了1e6个点来描述1~300区间内的z值,并将之存储为double双精度类型变量）；重新计算T1与TB；用T1-TB得出residu ;然后对于residu中凡是绝对值 < 容忍度(1e-5)的项目都认为是等于0.由此,我们便可以得出一系列在理论交点附近的点.最后再对这些点求平均值,就可以得出您的交点坐标了.

4.下面是我修改后的m文件.您直接复制粘贴后运行就行了.另外我额外做了一些图像处理并修改了一些您以前命令中的一些错误,希望能够帮助到您.

clear all

clc

% temperature diffusion 温度扩散曲线

z = 0:300;                                   %depth in mbsf 深度

t1 =  500;                                   % time

b = 0.054;                                   %thermal gradient

DT = 2;                                      %temperature change

k = 3E-7;                                    %thermal diffusivity

To = 2.2;                                      %initial temperature

T1 = To+b.*z+DT.*erfc(z./2.*sqrt(k.*t1)); % 温度随深度扩散曲线

%plot temperature diffusion line

figure('Name', 'Diagram for wulude','NumberTitle','off');

hold on

title('The intersection of 2 lines','FontSize',12);

plot(T1,z,'g','LineWidth',2);                               %green

xlabel('Temperature (C)','Color','b','FontWeight','bold');

ylabel('Depth (mbsf)','Color','b','FontWeight','bold');

grid on

% hydrate stability Line 稳定方程曲线

% 3.35% salt water 100% pure methane

P = 0.101325+1030.*9.81.*(750+z).*0.000001; %净水压力

LP = log10(P);

TB = 1./(0.00379-0.000283.*LP)-273.15;   %稳定方程曲线

% plot hydrate equilibrium

plot(TB,z,'m--','LineWidth',2);                       %megenta dash

%到这里都没问题

% junction point T z 求两条曲线交点,下面这行出错.

tolerance = 1e-5; 

% all the value which is smaller than this will be

% considered as 0

% here, in order to find a good result, we'll use a range more precise

% and recalculate the 2 functions

clear T1 TB z residu

z = linspace(1,300,1e6); % this time, we define 1e6 points between 1 and 300

T1 = To+b.*z+DT.*erfc(z./2.*sqrt(k.*t1));

P = 0.101325+1030.*9.81.*(750+z).*0.000001; %净水压力

LP = log10(P);

TB = 1./(0.00379-0.000283.*LP)-273.15;

residu = T1 - TB; % the difference between 2 functions

% some post-traitement of the result

% therefore we can locate the intersection point

results = [];

for i = 1 : length(z)

    if abs(residu(i)) < tolerance

        results = [results,i];

    end;

end;

l = length(results);

zResult = z(results(1):results(l));

tResult = T1(results(1):results(l));

zSol = mean(zResult);

tSol = mean(tResult);

% draw the point in the figure

plot(tSol,zSol,'r+','LineWidth',3,'MarkerSize',13);

disp('Reel Solution Found');

str = sprintf('The intersection of 2 lines \n z = %f   t = %f',zSol,tSol);

disp(str);

title(str,'FontSize',12);

 5.最后,由于我电脑装的是32位系统,所以我的最多只能用1e6个点来描述1~300这个区间.如果您使用的是64为系统+64位MATLAB,您可以增大这个数值,并且相应的减小容忍度,来无限接近理论交点——如果您的确需要一个十分十分精确的交点坐标

相反,您可以通过减少点数,增大容忍度,来快速定位您的交点——特别是当如果您处理多条有多个交点的函数时,此方法能够节省很多时间.并且满足您所要求的精确度.

希望上述两种求解方法都能帮助到您~

很高兴为您解答~