数学建模社区-数学中国

标题: 求助：关于分块矩阵的还原，为什么实行不了 [打印本页]
作者: 舒米牛牛    时间: 2009-12-17 20:17
标题: 求助：关于分块矩阵的还原，为什么实行不了
%运用Jacobi迭代拟合出来的u关于x和y的矩阵
0 G6 `/ e  @' Q; D3 T" ifunction J=jac(A,b,u0,eps)
0 U1 K: v+ J7 E9 H( Lif nargin==3
1 b6 d: z# ?! r( G5 [* r8 n1 H    eps=1.0e-8
' A: f; j7 P' |+ K- M  \: V- Gelseif nargin<3
: `; ^+ K8 ^( r    'error'; L. k8 V' E& Z# o
    return  U! a4 M/ ?6 X0 |' \6 s% Q
end  B/ a1 t( \4 n% l

/ p- S( V. o: u4 I: q8 e  B$ r%定义内部节点矩阵u09 w6 {" B% _& m8 C% }& `, p
h=1;k=1;' P' q0 J7 ^5 Y2 F6 ~, F
x=0:h:17;y=0:k:10;
2 J) l+ B; g( ?& h# Ye=length(x)-2;f=length(y)-2;8 H3 s- D; V! p) [% J0 ^* k0 ]
u=zeros(e,f);
% ]/ v; U  o# ^/ \$ }* T' ~u0=u;
" k7 J: S' r% B/ \5 b3 A7 Y6 T" y/ r4 r- a* k% b
%定义外部节点p6 \/ y; M$ _- U% n* c/ X+ l* |) ^
p=zeros(e+2,e+2);
! |) z; k" U# k5 m, Tp(1,1:f+2)=0;p(e+2,1:f+2)=0;       * Y* U, o+ D# D5 C0 W
p(1:e+2,1)=100;p(1:e+2,f+2)=100;
7 D; q4 s/ F4 T8 c0 q3 X) c$ N$ |
$ ]8 L3 ~' c: s/ A( c* O- Z%定义系数矩阵A- N2 r& |' L0 R$ y7 x' X3 C! w
A=zeros(e*f,e*f);' b+ Q8 b5 [) u5 J3 h5 U% r2 S
B=mat2cell(A,ones(e*f/e,1)*e,ones(e*f/e,1)*e);
' d; L- A2 G$ U) D8 S% ^d1=ones(e,1);d2=ones(e-1,1);
# b% }/ L8 s- p5 O: i. O/ b0 T. YM=4*diag(d1)-diag(d2,1)-diag(d2,-1);. e7 `) O: d& Q  Y( c
N=-eye(e);
$ K7 [2 I: D% [* JB{1,1}=M;B{e,e}=M
1 G; S  K' u6 Z9 w' t% _/ ]for i=2:e-1
/ \0 `4 [! P: ]% a    B{i,i}=M;B{i-1,i}=N;B{i+1,i}=N
/ [6 z+ m+ a- {- qend( H, \+ f* g) D) m
A=cell2mat(B);  
3 ]& J! W' S! ~! V- w0 n- I# B  m这里总是显示& g% \7 k5 x) k4 ^. d
??? function J=jac(A,b,u0,eps)
9 z" G3 u% \8 s. s2 ]7 f0 J    |
3 b/ C9 e0 L) Y0 LError: Function definitions are not permitted at the prompt or in scripts.
4 L: m0 B; Z2 k/ H- G. F  ?0 {. t/ c3 [
%定义b
" W4 Y, e8 ?# g( ]b=ones(e*f,1);
% u9 `' X7 S# r9 Tfor i=1:e5 W2 L5 d1 H% r  [) b
    for j=1:f
+ b; R. u7 W# V, E/ ]" Y4 {$ \) q! A        b(i+j)=p(i,j+1)+p(i+1,j)+p(i+2,j+1)+p(i+1,j+2)) `: F0 [( K, p, R% M
    end1 V/ a' e# Q7 M$ @' r1 U
end
2 P, S& C9 o/ O* X%运用Jacobi迭代法计算
! d5 W& I* I% l. N0 Y2 E. }D=diag(diag(A));2 t' E/ v& b0 l2 N8 S3 P' \) {, M
D=inv(D);; t7 c8 d/ p0 d$ v) z' O  T
L=tril(A,-1);
: I4 k( x( R  U; a) F8 wU=triu(A,1);
0 h- Y, Q# F9 c  q0 q7 h) D7 F8 uB=-D*(L+U);7 E* Y/ n: ]$ [1 q
f=D*b;
7 {1 i) N5 |" p3 DJ=B*u0+f;( Z# S' u) h/ ?% l) |% Y
while norm(J-u0)>=eps
2 R0 |* q/ v( W$ i3 q$ d  v4 K! B# Bx0=J;
- _% L( B5 v2 A, C3 P) y2 N" A! UJ=B*u0+f;
/ T. G- R! T2 }4 u+ u1 R; c# A. w" Xend* v8 b8 s) \# m' N- ^2 v
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