renminyingxong 发表于 2007-1-23 10:35

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

本人写了一个求解非线形方程的程序,系数可变,,有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 发表于 2007-1-23 11:09

使用reshape时候要注意变形前后的矩阵元素数量不能改变。
如X=,size(A)=,则X = reshape(X,size(A))之后有X=。
检查一下你的X和A吧。

renminyingxong 发表于 2007-1-23 14:14

回复 #2 realyyy 的帖子

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

realyyy 发表于 2007-1-24 10:40

你的pf和t1是什么?

renminyingxong 发表于 2007-1-24 11:30

回复 #4 realyyy 的帖子

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

realyyy 发表于 2007-1-24 13:26

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

xjzuo 发表于 2007-1-24 14:42

回复

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

renminyingxong 发表于 2007-1-24 15:40

回复 #7 xjzuo 的帖子

我的代码如下:
%一、将所用坐标值读入矩阵
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 发表于 2007-1-24 15:56

回复 #6 realyyy 的帖子

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

xjzuo 发表于 2007-1-24 16:25

回复

应该是清理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 发表于 2007-1-25 11:41

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

renminyingxong 发表于 2007-1-26 19:38

回复 #10 xjzuo 和#11 realyyy的帖子

感谢两位大虾的悉心讲解,我也将努力,此外,我还先就这个问题讨教一些其他的方面,希望两位大虾能够多多帮助.问题如下:
求解一个方程,其实就是求解一条直线与一条曲线的交点,原程序如下:
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
但是,我发现它给出的解很是奇怪,只给了一个解,而且很大,我将该直线与曲线在MATLAB中作图,发现它解不止一个啊,至少在250~300之间就有两个交点,即解.(不好意思,我画的图帖不上来).
我现在就是想请教一下,如何才能求出该方程的所有解,另外如果想求出一个范围内的解,应该怎么办,希望大虾能给个程序,我将详细拜读,再次致谢.

xjzuo 发表于 2007-1-26 21:53

回复

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

renminyingxong 发表于 2007-1-28 09:56

回复 #13 xjzuo 的帖子

谢谢您的回复,非常感谢
页: [1]
查看完整版本: 求助,望指点求解非线形方程的程序