数学建模社区-数学中国

标题: 求助：关于分块矩阵的还原，为什么实行不了 [打印本页]
作者: 舒米牛牛    时间: 2009-12-17 20:17
标题: 求助：关于分块矩阵的还原，为什么实行不了
%运用Jacobi迭代拟合出来的u关于x和y的矩阵
6 ]4 L9 k- U- afunction J=jac(A,b,u0,eps)
4 {1 i1 ~+ q( d! u1 A6 f4 E1 ]2 zif nargin==3( c/ J& ?' x0 E1 U' \* q
    eps=1.0e-8
9 m& P6 _  {: n2 N9 nelseif nargin<3
- D6 [1 i, k2 f  X6 \) u7 b& L    'error'% }$ S2 q1 {2 V( Z% y7 p! a. O
    return* ]( V9 A, O' d
end. p/ X7 \5 ~" b3 R

7 R, t! N( e0 w, h%定义内部节点矩阵u0
' R7 ?7 Q( Z  ^0 H8 ?h=1;k=1;
- I# ~  s3 z4 s: q% Ix=0:h:17;y=0:k:10;
6 z: E  F" C/ T1 u8 |8 Fe=length(x)-2;f=length(y)-2;" l" r: R- l. v  E- B7 R- ~
u=zeros(e,f);
/ ?, g- x4 O# y/ F# J/ }, {5 |u0=u;
8 F  l& a% K. e7 L4 M9 [7 i5 H" c
" |" ~& n7 ~. O6 }' @2 N+ A! M6 R%定义外部节点p4 z9 b9 T* o8 O4 B  u+ \+ P$ K
p=zeros(e+2,e+2);$ F3 T# j. g) e7 D4 W
p(1,1:f+2)=0;p(e+2,1:f+2)=0;       : D. e0 H) n# I( q7 f6 l
p(1:e+2,1)=100;p(1:e+2,f+2)=100;
5 g: L, @; V, g; s7 [0 I' M2 T4 o0 L# b5 M
%定义系数矩阵A
5 @& k7 \; p+ ^3 G& M: lA=zeros(e*f,e*f);- @2 W) v1 g/ Z* j2 |1 J/ v
B=mat2cell(A,ones(e*f/e,1)*e,ones(e*f/e,1)*e);
  Y+ ?5 F" z; ?- b8 _2 Ed1=ones(e,1);d2=ones(e-1,1);* h. i+ |0 B( U" h5 a5 |
M=4*diag(d1)-diag(d2,1)-diag(d2,-1);
& V% B. `2 f6 i/ v, ~N=-eye(e);  |* n3 E9 ?, m5 r' B
B{1,1}=M;B{e,e}=M
0 v) Y* Q. X5 ?+ p- \: c- T# H& _for i=2:e-18 P  {3 ~* z. Q7 ~% ~
    B{i,i}=M;B{i-1,i}=N;B{i+1,i}=N: d$ f' Z$ x. e3 J5 e3 b& r0 N- P
end; d5 B$ B% B3 G& m8 i
A=cell2mat(B);  ' l. v' m+ S! j' J* a
这里总是显示
1 P, @* a; s8 j4 O2 D??? function J=jac(A,b,u0,eps)* @' g8 w& i, W+ m9 z7 r
    |/ f. @- E' T6 M  o) _, z2 c  a
Error: Function definitions are not permitted at the prompt or in scripts.+ Q/ ~- G, j8 l7 p8 m
2 Z$ g/ T3 N, s# \: X0 s- H9 M! Y
%定义b7 _, s' u5 B2 Z7 M' {" n$ x# \
b=ones(e*f,1);6 }& _' `5 J4 V3 U# e! x* D
for i=1:e
/ n- U  `! {5 k- U    for j=1:f
4 G, E  f0 N7 I1 @- d) O. f! E        b(i+j)=p(i,j+1)+p(i+1,j)+p(i+2,j+1)+p(i+1,j+2)9 W& u4 r6 R+ P6 O! C$ V
    end) G- |3 o6 A- T8 j
end
/ D: E* O* G$ ]9 L3 Q( {) ~- H%运用Jacobi迭代法计算. P) ^& M' L. O' Z$ z( V& u# U
D=diag(diag(A));
' M2 a3 r" I; s( O; ]D=inv(D);
- s% E" T% ~" T  ?3 o. @% U- {) _L=tril(A,-1);
  l1 h7 q' c4 x! a# \+ T: a+ U$ J6 XU=triu(A,1);: @% |8 B- s  P  @& m$ B2 P( @
B=-D*(L+U);
. ^' |9 L6 ?* L1 W* @0 R) Ef=D*b;
3 Y: c$ |6 m1 W; }; I9 {: b. v( b! GJ=B*u0+f;
' d8 k: G" T" o/ Wwhile norm(J-u0)>=eps
) ^/ z4 I( ?* \1 I. ux0=J;( K( |) q' u; l' `1 {( e: q* F
J=B*u0+f;- m3 b- }" z. b, N
end0 b* @5 Z# q% w
return
作者: 舒米牛牛    时间: 2009-12-17 20:18
自己顶一下，拜托哪位高手指点一下，纠结这个矩阵的还原，想了好多方法还是不行~~实在想不出哪里出错了
作者: BenCam    时间: 2009-12-17 21:06
对不起,我也不知道,帮不了忙!
作者: madio    时间: 2009-12-17 22:38
你是不是函数的定义没有放在M文件中？



欢迎光临 数学建模社区-数学中国 (http://www.madio.net/) Powered by Discuz! X2.5