求助,望指点求解非线形方程的程序
本人写了一个求解非线形方程的程序,系数可变,,有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 编辑 ] 使用reshape时候要注意变形前后的矩阵元素数量不能改变。
如X=,size(A)=,则X = reshape(X,size(A))之后有X=。
检查一下你的X和A吧。
回复 #2 realyyy 的帖子
楼上的大虾,小弟没有看懂你说的意思,劳烦详细解释一下啊,我的程序中并没有使用RESHAPE啊,他内部调用这个出错,是怎么回事哈? 你的pf和t1是什么?回复 #4 realyyy 的帖子
我的pf和t1是两个数组,一维的,为这个问题我已经烦了好几天了,若大虾不烦,请加我QQ:66899343,淘金者.我将东西均发于你细看,还望不吝赐教 %在你的循环体内的q=solve(Y)的前面加入下面代码即可:Y=vpa(Y,4); %其中4可以根据需要调整;
if mod(n,500)==0 %其中500可以根据需要调整,越大越好。
clear maplemex
end
%速度飞快哦!
回复
建议LZ将问题和代码贴出来,一方面可以让大家也参与,另一方面也符合版规.回复 #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 编辑 ]
回复 #6 realyyy 的帖子
realyyy大师,我采用你的方法后,计算完成了,但只知其一不知其二,我在看到大师的帖子之前,采用了一个大循环,每循环1000次时就clear all一次,计算量很大,但也能计算下去,不如大师的算法精妙,在此很想知道大师的clear maplemex用意何在,望祥解,再次拜谢回复
应该是清理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 编辑 ] 我也是在大家的帮助下才知道的,clear maplemex是maple公司提供的一个折中解决matlab软件自身问题的命令,matlab在大量运算时有时会出现整数太大、内存不足等问题。
回复 #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之间就有两个交点,即解.(不好意思,我画的图帖不上来).
我现在就是想请教一下,如何才能求出该方程的所有解,另外如果想求出一个范围内的解,应该怎么办,希望大虾能给个程序,我将详细拜读,再次致谢.
回复
solve的确有时候不能给出全部解,所以建议用以下方式求解:%%%---------------------------------------------------------------%%%
对于你这种数值范围变化很大的问题,一般只能先画图,
再用fzero 或 fsolve求出相应的解.
%%%---------------------------------------------------------------%%%
回复 #13 xjzuo 的帖子
谢谢您的回复,非常感谢
页:
[1]