我用matlab编程计算扩散方程，答案总有负数，这个不正确，为什么
一，所要计算的扩散方程是：file:///C:/DOCUME~1/ADMINI~1/LOCALS~1/Temp/)EK`]U@Z(0~2D}6]_)$2SFA.jpg
这里c=0.164
a = 0.5
f = 1
边界条件， file:///C:/DOCUME~1/ADMINI~1/LOCALS~1/Temp/W01G9DX
g=0 ;q = 0.5;c = 1
file:///C:/DOCUME~1/ADMINI~1/LOCALS~1/Temp/~5%IIQ$D)K6)J`JI88X1LJF.jpg
file:///C:/DOCUME~1/ADMINI~1/LOCALS~1/Temp/54XFRZM_8TBXZR%~R2I1JYJ.jpg
file:///C:/DOCUME~1/ADMINI~1/LOCALS~1/Temp/S%@IOIWE04VJ03VIV0KTA9M.jpg
file:///C:/DOCUME~1/ADMINI~1/LOCALS~1/Temp/0KT@U7KMZ$TYVC1E6082_TB.jpg
file:///C:/DOCUME~1/ADMINI~1/LOCALS~1/Temp/QO$KQ14SHC2BDP)_N~~)4$B.jpg
file:///C:/DOCUME~1/ADMINI~1/LOCALS~1/Temp/$Y74]2TM0J.jpg
file:///C:/DOCUME~1/ADMINI~1/LOCALS~1/Temp/257POL2_DCILZ`]R9HOOWVA.jpg
file:///C:/DOCUME~1/ADMINI~1/LOCALS~1/Temp/FJHXZ)(C6BRY1ZXW)`124@B.jpg
file:///C:/DOCUME~1/ADMINI~1/LOCALS~1/Temp/T$9@BSHCX]9`SPWU2_83O75.jpg
二，问题是，所求出的u应该只有正值，现在却有很多负数，这个是不正确的。难道是我的代码编的有问题？？？？
三，代码实现如下
%初始化数组A，A的元素是点的坐标
A=double(zeros(331,2));
pre = 0;
for j=0:1:9%第几圈
x=0;
m=6*(j+1);
% fprintf(1,'m is %6.2f:\n',m);
for i=1:1:m%每一圈的点数
%fprintf(1,'x is %6.2f:\n',x);
k=i+pre;
A(k,1)=(j+1)*cos(x);
A(k,2)=(j+1)*sin(x);
x=x+pi/(3*(j+1));
% fprintf(1,'k is %6.2f:\n',k);
end
pre = pre+m;
%fprintf('\n');
end
%将圆心添加到数组里面
for i=331:-1:2
A(i,1)=A((i-1),1);
A(i,2)=A((i-1),2);
end
A(1,1)=0;A(1,2)=0;
%显示坐标的点
for i=1:1:331
fprintf(1,'第%6.2f个点：A(%6.2f,%6.2f)\n',i,A(i,1),A(i,2));
end
%画出所有的同心圆环
for i=1:1:10
t = -pi:0.01:pi;
x = i*cos(t);
y = i*sin(t);
plot(x,y);
hold on
end
%在圆上标出所取的点
for i=1:1:331
plot(A(i,1),A(i,2),'r*');
hold on
end
%定义B数组，这个数组的每个单元表示一个三角形，二维数组
B=[2,3,1;
2,9,3;
3,4,1;
3,11,4;
4,5,1;
4,13,5;
5,6,1;
5,15,6;
6,7,1;
6,17,7;
7,2,1;
7,19,2;
8,9,2;
8,21,9;
9,10,3;
9,22,10;
10,11,3;
10,24,11;
11,12,4;
11,25,12;
12,13,4;
12,27,13;
13,14,5;
13,28,14;
14,15,5;
14,30,15;
15,16,6;
15,31,16;
16,17,6;
16,33,17;
17,18,7;
17,34,18;
18,19,7;
18,36,19;
19,8,2;
19,37,8;
20,21,8;
20,39,21;
21,22,9;
21,40,22;
22,23,10;
22,41,23;
23,24,10;
23,43,24;
24,25,11;
24,44,25;
25,26,12;
25,45,26;
26,27,12;
26,47,27;
27,28,13;
27,48,28;
28,29,14;
28,49,29;
29,30,14;
29,51,30;
30,31,15;
30,52,31;
31,32,16;
31,53,32;
32,33,16;
32,55,33;
33,34,17;
33,56,34;
34,35,18;
34,57,35;
35,36,18;
35,59,36;
36,37,19;
36,60,37;
37,20,8;
37,61,20;
38,39,20;
38,63,39;
39,40,21;
39,64,40;
40,41,22;
40,65,41;
41,42,23;
41,66,42;
42,43,23;
42,68,43;
43,44,24;
43,69,44;
44,45,25;
44,70,45;
45,46,26;
45,71,46;
46,47,26;
46,73,47;
47,48,27;
47,74,48;
48,49,28;
48,75,49;
49,50,29;
49,76,50;
50,51,29;
50,78,51;
51,52,30;
51,79,52;
52,53,31;
52,80,53;
53,54,32;
53,81,54;
54,55,32;
54,83,55;
55,56,33;
55,84,56;
56,57,34;
56,85,57;
57,58,35;
57,86,58;
58,59,35;
58,88,59;
59,60,36;
59,89,60;
60,61,37;
60,90,61;
61,38,20;
61,91,38;
62,63,38;
62,93,63;
63,64,39;
63,94,64;
64,65,40;
64,95,65;
65,66,41;
65,96,66;
66,67,42;
66,97,67;
67,68,42;
67,99,68;
68,69,43;
68,100,69;
69,70,44;
69,101,70;
70,71,45;
70,102,71;
71,72,46;
71,103,72;
72,73,46;
72,105,73;
73,74,47;
73,106,74;
74,75,48;
74,107,75;
75,76,49;
75,108,76;
76,77,50;
76,109,77;
77,78,50;
77,111,78;
78,79,51;
78,112,79;
79,80,52;
79,113,80;
80,81,53;
80,114,81;
81,82,54;
81,115,82;
82,83,54;
82,117,83;
83,84,55;
83,118,84;
84,85,56;
84,119,85;
85,86,57;
85,120,86;
86,87,58;
86,121,87;
87,88,58;
87,123,88;
88,89,59;
88,124,89;
89,90,60;
89,125,90;
90,91,61;
90,126,91;
91,62,38;
91,127,62;
92,93,62;
92,129,93;
93,94,63;
93,130,94;
94,95,64;
94,131,95;
95,96,65;
95,132,96;
96,97,66;
96,133,97;
97,98,67;
97,134,98;
98,99,67;
98,136,99;
99,100,68;
99,137,100;
100,101,69;
100,138,101;
101,102,70;
101,139,102;
102,103,71;
102,140,103;
103,104,72;
103,141,104;
104,105,72;
104,143,105;
105,106,73;
105,144,106;
106,107,74;
106,145,107;
107,108,75;
107,146,108;
108,109,76;
108,147,109;
109,110,77;
109,148,110;
110,111,77;
110,150,111;
111,112,78;
111,151,112;
112,113,79;
112,152,113;
113,114,80;
113,153,114;
114,115,81;
114,154,115;
115,116,82;
115,155,116;
116,117,82;
116,157,117;
117,118,83;
117,158,118;
118,119,84;
118,159,119;
119,120,85;
119,160,120;
120,121,86;
120,161,121;
121,122,87;
121,162,122;
122,123,87;
122,164,123;
123,124,88;
123,165,124;
124,125,89;
124,166,125;
125,126,90;
125,167,126;
126,127,91;
126,168,127;
127,92,62;
127,169,92;
128,129,92;
128,171,129;
129,130,93;
129,172,130;
130,131,94;
130,173,131;
131,132,95;
131,174,132;
132,133,96;
132,175,133;
133,134,97;
133,176,134;
134,135,98;
134,177,135;
135,136,98;
135,179,136;
136,137,99;
136,180,137;
137,138,100;
137,181,138;
138,139,101;
138,182,139;
139,140,102;
139,183,140;
140,141,103;
140,184,141;
141,142,104;
141,185,142;
142,143,104;
142,187,143;
143,144,105;
143,188,144;
144,145,106;
144,189,145;
145,146,107;
145,190,146;
146,147,108;
146,191,147;
147,148,109;
147,192,148;
148,149,110;
148,193,149;
149,150,110;
149,195,150;
150,151,111;
150,196,151;
151,152,112;
151,197,152;
152,153,113;
152,198,153;
153,154,114;
153,199,154;
154,155,115;
154,200,155;
155,156,116;
155,201,156;
156,157,116;
156,203,157;
157,158,117;
157,204,158;
158,159,118;
158,205,159;
159,160,119;
159,206,160;
160,161,120;
160,207,161;
161,162,121;
161,208,162;
162,163,122;
162,209,163;
163,164,122;
163,211,164;
164,165,123;
164,212,165;
165,166,124;
165,213,166;
166,167,125;
166,214,167;
167,168,126;
167,215,168;
168,169,127;
168,216,169;
169,128,92;
169,217,128;
170,171,128;
170,219,171;
171,172,129;
171,220,172;
172,173,130;
172,221,173;
173,174,131;
173,222,174;
174,175,132;
174,223,175;
175,176,133;
175,224,176;
176,177,134;
176,225,177;
177,178,135;
177,226,178;
178,179,135;
178,228,179;
179,180,136;
179,229,180;
180,181,137;
180,230,181;
181,182,138;
181,231,182;
182,183,139;
182,232,183;
183,184,140;
183,233,184;
184,185,141;
184,234,185;
185,186,142;
185,235,186;
186,187,142;
186,237,187;
187,188,143;
187,238,188;
188,189,144;
188,239,189;
189,190,145;
189,240,190;
190,191,146;
190,241,191;
191,192,147;
191,242,192;
192,193,148;
192,243,193;
193,194,149;
193,244,194;
194,195,149;
194,246,195;
195,196,150;
195,247,196;
196,197,151;
196,248,197;
197,198,152;
197,249,198;
198,199,153;
198,250,199;
199,200,154;
199,251,200;
200,201,155;
200,252,201;
201,202,156;
201,253,202;
202,203,156;
202,255,203;
203,204,157;
203,256,204;
204,205,158;
204,257,205;
205,206,159;
205,258,206;
206,207,160;
206,259,207;
207,208,161;
207,260,208;
208,209,162;
208,261,209;
209,210,163;
209,262,210;
210,211,163;
210,264,211;
211,212,164;
211,265,212;
212,213,165;
212,266,213;
213,214,166;
213,267,214;
214,215,167;
214,268,215;
215,216,168;
215,269,216;
216,217,169;
216,270,217;
217,170,128;
217,271,170;
218,219,170;
218,273,219;
219,220,171;
219,274,220;
220,221,172;
220,275,221;
221,222,173;
221,276,222;
222,223,174;
222,277,223;
223,224,175;
223,278,224;
224,225,176;
224,279,225;
225,226,177;
225,280,226;
226,227,178;
226,281,227;
227,228,178;
227,283,228;
228,229,179;
228,284,229;
229,230,180;
229,285,230;
230,231,181;
230,286,231;
231,232,182;
231,287,232;
232,233,183;
232,288,233;
233,234,184;
233,289,234;
234,235,185;
234,290,235;
235,236,186;
235,291,236;
236,237,186;
236,293,237;
237,238,187;
237,294,238;
238,239,188;
238,295,239;
239,240,189;
239,296,240;
240,241,190;
240,297,241;
241,242,191;
241,298,242;
242,243,192;
242,299,243;
243,244,193;
243,300,244;
244,245,194;
244,301,245;
245,246,194;
245,303,246;
246,247,195;
246,304,247;
247,248,196;
247,305,248;
248,249,197;
248,306,249;
249,250,198;
249,307,250;
250,251,199;
250,308,251;
251,252,200;
251,309,252;
252,253,201;
252,310,253;
253,254,202;
253,311,254;
254,255,202;
254,313,255;
255,256,203;
255,314,256;
256,257,204;
256,315,257;
257,258,205;
257,316,258;
258,259,206;
258,317,259;
259,260,207;
259,318,260;
260,261,208;
260,319,261;
261,262,209;
261,320,262;
262,263,210;
262,321,263;
263,264,210;
263,323,264;
264,265,211;
264,324,265;
265,266,212;
265,325,266;
266,267,213;
266,326,267;
267,268,214;
267,327,268;
268,269,215;
268,328,269;
269,270,216;
269,329,270;
270,271,217;
270,330,271;
271,218,170;
271,331,218;
272,273,218;
273,274,219;
274,275,220;
275,276,221;
276,277,222;
277,278,223;
278,279,224;
279,280,225;
280,281,226;
281,282,227;
282,283,227;
283,284,228;
284,285,229;
285,286,230;
286,287,231;
287,288,232;
288,289,233;
289,290,234;
290,291,235;
291,292,236;
292,293,236;
293,294,237;
294,295,238;
295,296,239;
296,297,240;
297,298,241;
298,299,242;
299,300,243;
300,301,244;
301,302,245;
302,303,245;
303,304,246;
304,305,247;
305,306,248;
306,307,249;
307,308,250;
308,309,251;
309,310,252;
310,311,253;
311,312,254;
312,313,254;
313,314,255;
314,315,256;
315,316,257;
316,317,258;
317,318,259;
318,319,260;
319,320,261;
320,321,262;
321,322,263;
322,323,263;
323,324,264;
324,325,265;
325,326,266;
326,327,267;
327,328,268;
328,329,269;
329,330,270;
330,331,271;
331,272,218;]
%计算每个三角单元
L = cell(600,1);
a = double(zeros(331,331));%左边刚度矩阵1元素
d = double(zeros(331,1)); %右边一维向量
for i = 1:1:600
m = B(i,1);%i
n = B(i,2);%j
l = B(i,3);%k
b(m) = (A(n,2)-A(l,2)).*3./pi;%b（i）= div(Li,x) = (yj-yk)/2s
b(n) = (A(l,2)-A(m,2)).*3./pi;%b(j) = div(Lj,x) = (yk-yi)/2s
b(l) = (A(m,2)-A(n,2)).*3./pi;%b(k) = div(Lk,x) = (yi-yj)/2s
c(m) = (A(l,1)-A(n,1)).*3./pi;%c(i) = div(Li,y) = (xk-xj)/2s
c(n) = (A(m,1)-A(l,1)).*3./pi;%c(j) = div(Lj,y) = (xi-xk)/2s
c(l) = (A(n,1)-A(m,1)).*3./pi;%c(k) = div(Lk,y) = (xj-xi)/2s
%Li = [(xj.*yk-xk.*yj)+(yj-yk).*x+(xk-xj).*y]/2S
L(m,1) = {(A(n,1).*A(l,2)-A(l,1).*A(n,2)+(A(n,2)-A(l,2)).*x+(A(l,1)-A(n,1)).*y).*3/pi};
%Lj = [(xk.*yi-xi.*yk)+(yk-yi).*x+(xi-xk).*y]/2S
L(n,1) = {(A(l,1).*A(m,2)-A(m,1).*A(l,2)+(A(l,2))-A(m,2).*x+(A(m,1)-A(l,1)).*y).*3/pi};
%Lk = [(xi.*yj-xj.*yi)+(yi-yj).*x+(xj-xi).*y]/2S
L(l,1) = {(A(m,1).*A(n,2)-A(n,1).*A(m,2)+(A(m,2))-A(n,2).*x+(A(n,1)-A(m,1)).*y).*3./pi};
%三角形计算9次
for j = 1:1:3
M = B(i,j);%M是第i个三角单元，第j个元素
% fprintf(1,'M IS =====%6.2f\n',M);
for k = 1:1:3
N = B(i,k); %N是第i个三角单元，第k个元素
%二重积分计算
syms x y;
I = int(L{M,1}.*L{N,1},y,0,1-x);
P = int(I,x,0,1);
P = eval(P);%类型转换
%计算刚度矩阵
a(M,N) = a(M,N)+((b(M).*b(N)+c(M).*c(N)).*0.164./2+0.025.*P(1,1)).*pi./3;
end
end
% fprintf(1,'hahahahaha%6.2f\n',i);
end
%边界值计算
for i = 541:1:600
m = B(i,1);
n = B(i,2);
l = B(i,3);
b(m) = (A(n,2)-A(l,2)).*3./pi;%b（i）= div(Li,x) = (yj-yk)/2s
b(n) = (A(l,2)-A(m,2)).*3./pi;%b(j) = div(Lj,x) = (yk-yi)/2s
b(l) = (A(m,2)-A(n,2)).*3./pi;%b(k) = div(Lk,x) = (yi-yj)/2s
c(m) = (A(l,1)-A(n,1)).*3./pi;%c(i) = div(Li,y) = (xk-xj)/2s
c(n) = (A(m,1)-A(l,1)).*3./pi;%c(j) = div(Lj,y) = (xi-xk)/2s
c(l) = (A(n,1)-A(m,1)).*3./pi;%c(k) = div(Lk,y) = (xj-xi)/2s
%Li = [(xj.*yk-xk.*yj)+(yj-yk).*x+(xk-xj).*y]/2S
L(m,1) = {(A(n,1).*A(l,2)-A(l,1).*A(n,2)+(A(n,2)-A(l,2)).*x+(A(l,1)-A(n,1)).*y).*3./pi};
%Lj = [(xk.*yi-xi.*yk)+(yk-yi).*x+(xi-xk).*y]/2S
L(n,1) = {(A(l,1).*A(m,2)-A(m,1).*A(l,2)+(A(l,2))-A(m,2).*x+(A(m,1)-A(l,1)).*y).*3./pi};
%Lk = [(xi.*yj-xj.*yi)+(yi-yj).*x+(xj-xi).*y]/2S
L(l,1) = {(A(m,1).*A(n,2)-A(n,1).*A(m,2)+(A(m,2))-A(n,2).*x+(A(n,1)-A(m,1)).*y).*3./pi};
for j = 1:1:3
M = B(i,j);%M是第i个三角单元，第j个元素
if j<=2
syms x y;
y = 0;
I2 = int(L{M,1}.*L{N,1},x,0,1);
P2 = eval(I2);
a(M,N) = a(M,N)+pi./6.*P2(1,1);
% fprintf(1,'a(%d,%d)=%6.2f\t\t\t\t',M,N,a(M,N));
end
%计算光源的那一项
if j==2&mod(i,5)==0
syms x y;
I1 = int(3.*L{M,1},y,0,1-x);
P1 = int(I1,x,0,1);
P1 = eval(P1);
%fprintf(1,'P1 =========is %6.2f\n',P1(1,1));
%fprintf(1,'MMMMMMM======%6.2f\n',M(1,1));
d(M) =d(M)+P1(1,1).*pi./3;
end
end
end
u= inv(a)*d;
呵呵，现在终于有答案了，主要原因是求a时出了些小错误 呵呵，幸好自己找出错误了，不然这么长的程序，想差错可费劲啊！
建议分享一下自己的处理过程！ 没错, 这么长的程序, 想查错可能很费劲啊!
建议LZ下次简化下, 并看下以下连接:@)
建议提问的网友分清 编程问题 和 专业问题
http://forum.vibunion.com/forum/viewthread.php?tid=36746&extra=&page=1
提问的智慧！！！！(发帖前请认真阅读)
http://forum.vibunion.com/forum/viewthread.php?tid=21991 回复 ChaChing 的帖子
这个问题有点困难，但是，是其中的一项出错了。 检查程序错误是一个头疼的问题，记得我们老师有一次在程序中+号写成-号，检查一个月才查出来，我都是打印出来一步步的查的 源码传下
页:
[1]