新手向您请教一简单问题，谢谢！

airhoo1

文本文档1.txt为（7列12行）：
0.19 0.07 0.74 0.600957 1.5593 0.4159 13.18293872
0.19 0.07 0.74 0.61483   1.5644 0.3932 13.18977909
0.19 0.07 0.74 0.563157 1.5672 0.3764 13.18983356
0.19 0.07 0.74 0.580532 1.5627 0.3819 13.1967284
0.19 0.07 0.74 0.572817 1.5485 0.385 13.19683803
0.19 0.07 0.74 0.574459 1.5437 0.3827 13.19683803
0.19 0.07 0.74 0.558511 1.5324 0.4194 13.21719435
0.19 0.07 0.74 0.637349 1.5381 0.4576 13.21719435
0.19 0.07 0.74 0.617176 1.5593 0.4545 13.21057324
0.19 0.07 0.74 0.634689 1.5526 0.4435 13.21057324
0.19 0.07 0.74 0.532446 1.5526 0.4141 13.23120398
0.19 0.07 0.74 0.517933 1.5567 0.3781 13.23120398

我要将文本文件1.txt每行的数值赋值给7个变量，将这些变量带入方程后，解关于Sh的一个一元多次方程，最后返回计算出的Sh结果，应该得到一个一维列向量。我试着将第一行的数值赋值给7个变量后计算Sh，但由于初学，总调试不正确。由于工作紧迫，特向您请教我的错误该如何更正？另外，要使每行都计算出一个Sh，是否将程序中第二行变为i=1:12？感谢您的赐教！谢谢！

我写的程序为：
f=load('1.txt');
i=1;
c1=f(i,1);c2=f(i,2);c3=f(i,3);o=f(i,4);Vp=f(i,5);Vs=f(i,6);p=f(i,7);
Km=0.5*(c1*36.6+c2*76+c3*6.85+(c1/36.6+c2/76+c3/6.85)^(-1));
Um=0.5*(c1*45+c2*26+c3*20.9+(c1/45+c2/26+c3/20.9)^(-1));
a=sqrt(Um/(Km+4/3*Um));
G=0.9552+0.0448*exp(-c3/0.06714);
m=1;
n=10^(0.426-0.235*log10(p))/m;
sym Sh;
y=Vp*(G*a*(1-o)^n+sqrt(3.3*Sh/(((1-Sh)/2.29+Sh/7.9)^(-1)+4/3*3.3*Sh))*o*Sh)-Vs;
Sh=solve(y,Sh)

显示错误是：
??? Error using ==> sym.mpower
Matrix must be square.

我用的matlab 7.0版本。

感激不尽！

sigma665

你最好逐行运行，这样才知道哪一行出问题，目前看来是乘有问题，不是方阵
另外，请搜索下精华帖关于文件读取，里面有很详细说明

airhoo1

谢谢sigma665老师，我逐行试了一下，是在倒数第二行时出错。请教一下该如何改正？谢谢！

ChaChing

...
显示错误是：
??? Error using ==> sym.mpower
Matrix must be square. ...

我报错怎不同, 楼主是否搞乱了!
??? Undefined function or variable 'Sh'.

试试将sym Sh改为syms Sh

airhoo1

sorry，是搞乱了，报错为??? Undefined function or variable 'Sh'

我将sym改为syms Sh后，运行得到的是一个表达式，而不是一个数值啊，请问是怎么回事？该怎么改正呢？

谢谢！

ChaChing

4F不是给了建议了!

airhoo1

谢谢各位的回答，我在最后加了一行subs(Sh)后已能计算出数值
f=load('1.txt');
i=1;
c1=f(i,1);c2=f(i,2);c3=f(i,3);o=f(i,4);Vp=f(i,5);Vs=f(i,6);p=f(i,7);
Km=0.5*(c1*36.6+c2*76+c3*6.85+(c1/36.6+c2/76+c3/6.85)^(-1));
Um=0.5*(c1*45+c2*26+c3*20.9+(c1/45+c2/26+c3/20.9)^(-1));
a=sqrt(Um/(Km+4/3*Um));
G=0.9552+0.0448*exp(-c3/0.06714);
m=1;
n=10^(0.426-0.235*log10(p))/m;
Sh=sym('Sh');
y=Vp*(G*a*(1-o)^n+sqrt(3.3*Sh/(((1-Sh)/2.29+Sh/7.9)^(-1)+4/3*3.3*Sh))*o*Sh)-Vs;
Sh=solve(y,Sh);
subs(Sh)

现在我想计算1.txt中1到12行所有的Sh，将第2行改为i=1:12；后出错??? Error using ==> mpower
Matrix must be square.，请问如何实现这个循环呢？谢谢！

ChaChing

f=load('1.txt');
for i=1:12
c1=f(i,1);c2=f(i,2);c3=f(i,3);o=f(i,4);Vp=f(i,5);Vs=f(i,6);p=f(i,7);
Km=0.5*(c1*36.6+c2*76+c3*6.85+(c1/36.6+c2/76+c3/6.85)^(-1));
Um=0.5*(c1*45+c2*26+c3*20.9+(c1/45+c2/26+c3/20.9)^(-1));
a=sqrt(Um/(Km+4/3*Um));
G=0.9552+0.0448*exp(-c3/0.06714);
m=1;
n=10^(0.426-0.235*log10(p))/m;
syms Sh;
y=Vp*(G*a*(1-o)^n+sqrt(3.3*Sh/(((1-Sh)/2.29+Sh/7.9)^(-1)+4/3*3.3*Sh))*o*Sh)-Vs;
Sh=solve(y,Sh);
subs(Sh)
end

[ 本帖最后由 ChaChing 于 2009-2-12 13:24 编辑 ]

airhoo1

非常感谢！得出了计算结果
ans =
    1.4014
    0.2876
ans =
    0.2657
    1.4027
ans =
    0.1404
    1.4072 ...

12行，每行得出Sh有2个解。我还想请教一下如何使计算结果显示为一个2列12行的数据表，因为我真正计算时有一千多行，如果人工整理上面这种形式的计算结果太费时间了。麻烦了！谢谢！

[ 本帖最后由 ChaChing 于 2009-2-12 20:37 编辑 ]

airhoo1

我的意思是使输出的结果变成两列，例如上面的结果变为下列形式：
ans=
1.4014 0.2876
1.4027 0.2657
1.4072 0.1404
1.4059 0.1890
1.4060 0.1876
1.4060 0.1885
1.4040 0.2380
1.3908 0.4146
1.3947 0.3757
1.3941 0.3814

谢谢您的帮助！

ChaChing

n=size(f,1); TolSh=zeros(n,2);
for i=1:n
c1=f(i,1);c2=f(i,2);c3=f(i,3);o=f(i,4);Vp=f(i,5);Vs=f(i,6);p=f(i,7); Km=0.5*(c1*36.6+c2*76+c3*6.85+(c1/36.6+c2/76+c3/6.85)^(-1));
Um=0.5*(c1*45+c2*26+c3*20.9+(c1/45+c2/26+c3/20.9)^(-1)); a=sqrt(Um/(Km+4/3*Um)); G=0.9552+0.0448*exp(-c3/0.06714);
m=1; n=10^(0.426-0.235*log10(p))/m; syms Sh; y=Vp*(G*a*(1-o)^n+sqrt(3.3*Sh/(((1-Sh)/2.29+Sh/7.9)^(-1)+4/3*3.3*Sh))*o*Sh)-Vs;
Sh=solve(y,Sh); TolSh(i,:)=subs(Sh)';
end
TolSh

airhoo1

谢谢ChaChing！可能有一些无解，运行后出现错误：
Warning: Explicit solution could not be found.
> In solve at 140
  In sym.solve at 49
??? Function 'eval' is not defined for values of class 'double'.

Error in ==> eval at 44
  [varargout{1:nargout}] = builtin('eval', varargin{:});

Error in ==> sym.double at 45
D = reshape(eval(X),siz);

Error in ==> sym.subs at 66
if isempty(vars), NEWf = double(OLDf); return, end
请问该如何解决呢？谢谢！

ChaChing

symbolic这块个人以前根本不会!
我是执行网友的程序, 边看Help, 边学习的!
LZ这样的讯息我就没经验了! 除非我可以试!
或等有经验的人出现! sorry!

刚刚才看到LZ数据! 有空会试试看!

[ 本帖最后由 ChaChing 于 2009-2-13 10:36 编辑 ]

airhoo1

:lol 还是感谢你！其实我也应该像你边执行边看help学习，可惜是初接触matlab，没基础，加上现在急需得出结果，所以一步步程序的实现都麻烦大家了。我自己也在摸索着，以后还要向大家多多学习！再次感谢你，并期待着能解决这个问题的人出现，感谢大家的帮助！

airhoo1

另外，我把程序中的文件1.txt发短消息给你了，也期待着你能攻克这个问题，我现在也在试着。麻烦你了！Thanks！