求助,望指点求解非线形方程的程序

renminyingxong

本人写了一个求解非线形方程的程序,系数可变,,有430000个数据,所以方程需要循环来求解,,程序如下:
syms x
for n=1:number
    x1=273.15;y1=pf(n);x2=t1(n);y2=0.375;
    a=(y2-y1)/(x2-x1);b=y1-a*x1;
    Y=1000*(a*x+b)-exp(14.717-1886.79/x);
    q=solve(Y);t=double(q);p=a*t+b;
end
number有430000个数,但是他在循环计算了3000多次后,就报错了,如下:
??? Error using ==> reshape
To RESHAPE the number of elements must not change.

Error in ==> sym.minus at 29
X = reshape(X,size(A));

Error in ==> shuzhijisuan2 at 27
    p1=1000*(a*x+b)-exp(14.717-1886.79/x);
请高手指点以下哈,小弟将不胜感激

[ 本帖最后由 mjhzhjg 于 2007-1-24 23:48 编辑 ]

realyyy

使用reshape时候要注意变形前后的矩阵元素数量不能改变。
如X＝[1,2,3,4,5,6]，size(A)＝[2,3]，则X = reshape(X,size(A))之后有X＝[1,3,5;2,4,6]。
检查一下你的X和A吧。

renminyingxong

楼上的大虾,小弟没有看懂你说的意思,劳烦详细解释一下啊,我的程序中并没有使用RESHAPE啊,他内部调用这个出错,是怎么回事哈?

realyyy

你的pf和t1是什么？

renminyingxong

我的pf和t1是两个数组,一维的,为这个问题我已经烦了好几天了,若大虾不烦,请加我QQ:66899343,淘金者.我将东西均发于你细看,还望不吝赐教

realyyy

%在你的循环体内的q=solve(Y)的前面加入下面代码即可：
    Y=vpa(Y,4);   %其中4可以根据需要调整；
    if mod(n,500)==0   %其中500可以根据需要调整，越大越好。
        clear maplemex
    end
%速度飞快哦!

xjzuo

建议LZ将问题和代码贴出来，一方面可以让大家也参与，另一方面也符合版规．

renminyingxong

我的代码如下:
%一、将所用坐标值读入矩阵
number=437437;%经纬度的个数，也就是所用文本文件的行数
jingdu=zeros(1,number);weidu=zeros(1,number);
pf=zeros(1,number);ps=zeros(1,number);t1=zeros(1,number);
guodu=zeros(3,1);
fid=fopen('d:\pf.txt','r');
for n=1:number
    guodu=fscanf(fid,'%f %f %f\n',3);jingdu(n)=guodu(1);
    weidu(n)=guodu(2);pf(n)=guodu(3);
end
fclose(fid);
fid=fopen('d:\t1.txt','r');
for n=1:number
    guodu=fscanf(fid,'%f %f %f\n',3);t1(n)=guodu(3);
end
fclose(fid);
%二、计算AB直线与曲线的交点，并写入AB文本
syms x
fid=fopen('d:\AB.txt','w');
for n=1:number
    x1=273.15;y1=pf(n);x2=t1(n);y2=0.375;
    a=(y2-y1)/(x2-x1);b=y1-a*x1;
    Y=1000*(a*x+b)-exp(14.717-1886.79/x);
    q=solve(Y);t=double(q);p=a*t+b;
    if imag(t)~=0
       t=0;
    end
    if imag(p)~=0
       p=0;
    end
    fprintf(fid,'%f  %f  %f\n',jingdu(n),weidu(n),t);
    fprintf(fid,'%f  %f  %f\n',jingdu(n),weidu(n),p);
end
fclose(fid);
pf文本很大,共有437437行,部分如下:
97.2866        39.7036        3.16
97.295        39.7036        2.74
97.3033        39.7036        2.66
97.32        39.6953        2.74
97.37        39.6786        2.44
97.37        39.6703        2.56
97.3783        39.6703        2.42
97.42        39.6703        2.49
97.4283        39.6703        2.96
97.4366        39.6703        2.63
97.395        39.662        2.50
97.4033        39.662        2.43
97.4116        39.662        2.64
96.845        39.6453        2.52
96.8533        39.6453        2.64
97.57        39.637        2.68
97.5283        39.6286        2.46
97.5616        39.6286        2.44
97.57        39.6286        2.71
97.5783        39.6286        2.69
97.5866        39.6286        2.74
97.52        39.6203        2.76
97.5283        39.6203        2.94
97.5366        39.6203        2.45
96.3033        39.612        2.64
96.345        39.6036        2.67
96.3533        39.6036        2.52
96.3367        39.5953        2.50
96.345        39.5953        2.44
96.3533        39.5953        2.76
95.2451        39.587        2.48
95.2534        39.587        2.60
95.3367        39.587        2.58
95.345        39.587        2.57
95.4284        39.587        2.47
96.22        39.587        2.52
96.2283        39.587        2.86
96.2367        39.587        3.31
96.245        39.587        3.39
96.2533        39.587        2.43
96.27        39.587        2.53
96.2783        39.587        2.61
96.2867        39.587        2.90
96.3033        39.587        2.60
96.3117        39.587        2.99
96.32        39.587        3.10
96.3283        39.587        3.33
96.3367        39.587        3.45
96.345        39.587        2.98
96.3533        39.587        3.27
96.3617        39.587        2.52
96.4033        39.587        2.39
95.2451        39.5786        2.80
95.2534        39.5786        3.11
95.2617        39.5786        2.69
t1文本也很大,共有437437行,部分如下:
97.2866        39.7036        268.44
97.295        39.7036        269.15
97.3033        39.7036        269.27
97.32        39.6953        269.14
97.37        39.6786        269.65
97.37        39.6703        269.45
97.3783        39.6703        269.68
97.42        39.6703        269.57
97.4283        39.6703        268.77
97.4366        39.6703        269.32
97.395        39.662        269.55
97.4033        39.662        269.66
97.4116        39.662        269.32
96.845        39.6453        269.52
96.8533        39.6453        269.30
97.57        39.637        269.24
97.5283        39.6286        269.62
97.5616        39.6286        269.64
97.57        39.6286        269.20
97.5783        39.6286        269.23
97.5866        39.6286        269.15
97.52        39.6203        269.12
97.5283        39.6203        268.82
97.5366        39.6203        269.64
96.3033        39.612        269.30
96.345        39.6036        269.27
96.3533        39.6036        269.52
96.3367        39.5953        269.54
96.345        39.5953        269.64
96.3533        39.5953        269.10
95.2451        39.587        269.58
95.2534        39.587        269.39
95.3367        39.587        269.41
95.345        39.587        269.44
95.4284        39.587        269.59
96.22        39.587        269.52
96.2283        39.587        268.93
96.2367        39.587        268.19
96.245        39.587        268.05
96.2533        39.587        269.66
96.27        39.587        269.49
96.2783        39.587        269.37
96.2867        39.587        268.87
96.3033        39.587        269.37
96.3117        39.587        268.72
96.32        39.587        268.53
96.3283        39.587        268.15
96.3367        39.587        267.95
96.345        39.587        268.74
96.3533        39.587        268.26
96.3617        39.587        269.51
96.4033        39.587        269.73
95.2451        39.5786        269.05
95.2534        39.5786        268.52
95.2617        39.5786        269.23
95.27        39.5786        269.53
95.2784        39.5786        269.42
95.2867        39.5786        269.57
95.345        39.5786        269.18
95.3534        39.5786        269.28
95.3617        39.5786        269.10
95.3784        39.5786        269.39
95.6534        39.5786        269.42
95.6617        39.5786        269.22
现在的问题是:计算到1000多行时就报错,而且是随即的,我试算多次,有时算到1200多行,有时算到1400多行,总之就是不确定.不知道问题出在哪里.
报错如下:
??? Error using ==> reshape
To RESHAPE the number of elements must not change.

Error in ==> sym.minus at 29
X = reshape(X,size(A));

Error in ==> shuzhijisuan2 at 27
    p1=1000*(a*x+b)-exp(14.717-1886.79/x);
这个就是整个问题,感谢斑竹的提醒,在此道歉先

[ 本帖最后由 renminyingxong 于 2007-1-24 15:48 编辑 ]

renminyingxong

realyyy大师,我采用你的方法后,计算完成了,但只知其一不知其二,我在看到大师的帖子之前,采用了一个大循环,每循环1000次时就clear all一次,计算量很大,但也能计算下去,不如大师的算法精妙,在此很想知道大师的clear maplemex用意何在,望祥解,再次拜谢

xjzuo

应该是清理Maple的工作空间.
目的是避免每次运算时都执行一些重复的过程．
建议看看http://forum.vibunion.com/forum/viewthread.php?tid=23984&highlight=%BC%C6%CB%E3%B9%FD%B3%CC%D6%D0%D5%FB%CA%FD%CC%AB%B4%F3

另:不少用Matlab7.0以上版本的人似乎也会常碰到此类错误提示,
     希望有详细研究过此类问题的人可以作一次较详细地分析和讲解.

[ 本帖最后由 xjzuo 于 2007-1-24 17:56 编辑 ]

realyyy

我也是在大家的帮助下才知道的，clear maplemex是maple公司提供的一个折中解决matlab软件自身问题的命令，matlab在大量运算时有时会出现整数太大、内存不足等问题。

renminyingxong

感谢两位大虾的悉心讲解，我也将努力，此外，我还先就这个问题讨教一些其他的方面，希望两位大虾能够多多帮助．问题如下：
求解一个方程，其实就是求解一条直线与一条曲线的交点，原程序如下：
syms x
x1=273.15;y1=3.16;x2=275.15;y2=4.14;
a=(y2-y1)/(x2-x1);b=y1-a*x1;
Y=1000*(a*x+b)-exp(38.98-8533.8/x);
S=solve(Y)
该程序运行结果如下：
>> syms x
x1=273.15;y1=3.16;x2=275.15;y2=4.14;
a=(y2-y1)/(x2-x1);b=y1-a*x1;
Y=1000*(a*x+b)-exp(38.98-8533.8/x);
S=solve(Y)

S =

173221910629935.53947320850193681
但是，我发现它给出的解很是奇怪，只给了一个解，而且很大，我将该直线与曲线在ＭＡＴＬＡＢ中作图，发现它解不止一个啊，至少在２５０～３００之间就有两个交点，即解．（不好意思，我画的图帖不上来）．
我现在就是想请教一下，如何才能求出该方程的所有解，另外如果想求出一个范围内的解，应该怎么办，希望大虾能给个程序，我将详细拜读，再次致谢．

xjzuo

solve的确有时候不能给出全部解，所以建议用以下方式求解：
%%%---------------------------------------------------------------%%%
对于你这种数值范围变化很大的问题,一般只能先画图,
再用fzero 或 fsolve求出相应的解．
%%%---------------------------------------------------------------%%%

renminyingxong

谢谢您的回复,非常感谢