数学建模社区-数学中国

标题: 求助：关于分块矩阵的还原，为什么实行不了 [打印本页]
作者: 舒米牛牛    时间: 2009-12-17 20:17
标题: 求助：关于分块矩阵的还原，为什么实行不了
%运用Jacobi迭代拟合出来的u关于x和y的矩阵
" Z+ o+ x2 l6 _4 z5 Bfunction J=jac(A,b,u0,eps)
3 u! z$ b& h! M- O8 Tif nargin==3; V* I7 I3 w* a
    eps=1.0e-8( e3 C3 `, A9 i; H& h3 B6 m
elseif nargin<3
/ Q' X( Z7 t/ d+ S0 C* j7 }    'error'; ^9 A. r5 ~: w' f0 J0 p
    return; O7 x9 ]! t; B* T/ `6 {+ t
end% S1 R7 X( `7 L& V6 |3 ?1 ]
- z9 I( J- B& @+ y2 m+ ?! R
%定义内部节点矩阵u05 d' F5 `! n7 K  O, C
h=1;k=1;
% N0 k8 ?. F3 h# s3 `2 rx=0:h:17;y=0:k:10;7 y7 q8 B: j' Q
e=length(x)-2;f=length(y)-2;; H- k" q3 ^8 F) L7 ^7 Y
u=zeros(e,f);
! B& }3 O: H' [8 E0 o: Tu0=u;
  L8 V: @' I7 E# @3 u/ o& t
, w2 S& a" b  ]7 E; w5 P: @%定义外部节点p+ W* r. T: i$ ?6 J
p=zeros(e+2,e+2);
7 w# z6 k( q; \p(1,1:f+2)=0;p(e+2,1:f+2)=0;         |7 c2 ]6 K: _  c  \3 ?
p(1:e+2,1)=100;p(1:e+2,f+2)=100;8 C. K. `1 C& g( h# Y. R

% k8 p1 M: k1 W( M%定义系数矩阵A
' F3 J, ?' i1 ^! Y! w5 I; HA=zeros(e*f,e*f);7 j, f, G. ]  j. N% Y2 P
B=mat2cell(A,ones(e*f/e,1)*e,ones(e*f/e,1)*e);
* w. Q) N! ?  Y' h3 }d1=ones(e,1);d2=ones(e-1,1);/ B8 T! T* _) X
M=4*diag(d1)-diag(d2,1)-diag(d2,-1);
2 F! z9 _% i$ @+ |  M1 gN=-eye(e);; i" D- ~9 d% ]% J& ~3 D: C
B{1,1}=M;B{e,e}=M
; \7 I# i7 k, E0 r& t' c* y7 }5 D+ dfor i=2:e-17 A  X0 i4 a5 Z1 E' Z; U8 ]
    B{i,i}=M;B{i-1,i}=N;B{i+1,i}=N
3 a$ O# k! w, B# J/ a. ~# }end2 _+ P7 N* P- }; y6 n6 X
A=cell2mat(B);  
8 E  }0 f' R# c这里总是显示# z9 X' q0 Z& Z
??? function J=jac(A,b,u0,eps)# {7 x) _9 o: V
    |/ x2 G3 g" t  p2 ?
Error: Function definitions are not permitted at the prompt or in scripts.
% d/ Y( |. {- @1 O& s* r" P& [* N8 }, f/ Q* L
%定义b
; m' H* r7 E; z) \, B; ^b=ones(e*f,1);( s5 w5 G5 _5 z2 l6 b' t! R& I
for i=1:e
( @$ i5 R0 W- h    for j=1:f! d# @3 i8 e( ?4 E
        b(i+j)=p(i,j+1)+p(i+1,j)+p(i+2,j+1)+p(i+1,j+2)
! o( Q( g  ]# O! M' h    end" w1 p, x3 r4 d3 [  y0 B
end5 z+ N7 Y1 {/ H
%运用Jacobi迭代法计算, p9 i( W- C7 |$ d
D=diag(diag(A));
6 {. z- k5 I; ^) J( nD=inv(D);
+ x+ j/ }' Q; B3 rL=tril(A,-1);) w- O" s: A, a, q$ W! b  Y
U=triu(A,1);
1 V% q) x3 R. g9 D5 ZB=-D*(L+U);" j6 C( R& {2 J. I8 {- m% i) D7 a
f=D*b;8 o7 U& i/ y/ a- B( G, z7 j/ N& K6 e
J=B*u0+f;
) ^5 y1 j0 {5 K4 Cwhile norm(J-u0)>=eps
& E0 k* q" [' L$ O, qx0=J;
  }8 X- ?* Q+ J1 G3 F2 _J=B*u0+f;" D; i" t7 E' K" M* m4 c
end
5 X1 z1 m- p  L: jreturn
作者: 舒米牛牛    时间: 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