数学建模社区-数学中国

标题: 求助：关于分块矩阵的还原，为什么实行不了 [打印本页]
作者: 舒米牛牛    时间: 2009-12-17 20:17
标题: 求助：关于分块矩阵的还原，为什么实行不了
%运用Jacobi迭代拟合出来的u关于x和y的矩阵
- B. R7 e$ p" e: C0 Ffunction J=jac(A,b,u0,eps)
; O6 [; k( T$ r' j  r+ P, xif nargin==30 K5 N- X8 X' {7 O
    eps=1.0e-8
0 Y" j, r- V( R9 n% h+ k+ Zelseif nargin<31 m# Q) f% u0 @) u7 h: {
    'error'
4 v- q( R0 W6 |3 `5 N    return
/ a) C% ]/ N! Q2 H& Send+ s, J- a* X+ D

0 W' q! ^; s, H7 A3 A6 Q* e4 h%定义内部节点矩阵u0( U1 m/ x1 z4 x5 f
h=1;k=1;
' ?& x2 m; I7 m6 |0 }+ n2 W- i8 Ax=0:h:17;y=0:k:10;# U2 v$ m( G7 E, d3 G
e=length(x)-2;f=length(y)-2;
2 p% J1 y3 }* R% O, N. N' i; X' l5 Ru=zeros(e,f);
7 J; p$ G9 i- d6 fu0=u;9 C$ ?$ g" ?- R& W+ t1 W) G
: a1 R" L+ K# R2 U+ ]7 A
%定义外部节点p  y  t4 X7 b) N
p=zeros(e+2,e+2);
7 ^5 i, @# _& Gp(1,1:f+2)=0;p(e+2,1:f+2)=0;      
1 d# T; v( B: L8 _# v0 L: c& Wp(1:e+2,1)=100;p(1:e+2,f+2)=100;
# J  R  S/ I4 \1 z+ n" |2 X
1 q. L6 C& |  b2 Z%定义系数矩阵A8 D3 j. a* N& I
A=zeros(e*f,e*f);: A. j2 P( T( r, @$ f0 h0 d
B=mat2cell(A,ones(e*f/e,1)*e,ones(e*f/e,1)*e);% A2 N0 R9 I" f: f; S- {. d% k
d1=ones(e,1);d2=ones(e-1,1);* O7 r5 z5 h% [0 A) h7 u
M=4*diag(d1)-diag(d2,1)-diag(d2,-1);9 Z) u0 X- U9 @  L4 K+ @
N=-eye(e);
) @+ }( ?( }$ T6 y1 t& s3 y, SB{1,1}=M;B{e,e}=M
+ G& G+ E; @8 X) v2 o# _7 dfor i=2:e-1( E" A. N3 _! X" `& _) ]. [
    B{i,i}=M;B{i-1,i}=N;B{i+1,i}=N
9 l4 R0 R. o) c# b5 Z1 z# L7 r( [end6 ?; U& e; ?/ i6 r3 p
A=cell2mat(B);  
6 R$ H- X5 W2 s2 d5 B0 ^. L  k9 k" F1 d4 T这里总是显示9 H  }1 M4 O& J
??? function J=jac(A,b,u0,eps)
5 K+ s( v7 }3 A! Q) X0 p2 ]    |8 C% H: w% H/ ^/ k" _0 i
Error: Function definitions are not permitted at the prompt or in scripts.
2 u/ S1 j& E9 l" R( M
, v& O9 P& b4 h' a/ T- I% F%定义b0 u7 N" [( h: J" I, G* _3 |. N! h
b=ones(e*f,1);/ j( f. N1 m) \& X% u% a% o# Q
for i=1:e1 V; \. f# m, z  z
    for j=1:f
6 y7 ~/ Q( ?+ b1 n        b(i+j)=p(i,j+1)+p(i+1,j)+p(i+2,j+1)+p(i+1,j+2)
; q3 ~) l' X' m    end
4 ~) R4 ^# `, G6 l3 E( U. r% ~end3 e/ i3 w+ t, @. w6 M
%运用Jacobi迭代法计算
4 i6 |% `' a/ K8 G; S- [* x. R; lD=diag(diag(A));
7 t' }3 y4 c0 P0 O0 Y: x/ s0 ID=inv(D);
; w, [6 c5 K/ ~L=tril(A,-1);
3 O% L: t) ?& OU=triu(A,1);( u. w9 U' L1 J. U* C
B=-D*(L+U);
: J! Y% J- [1 [f=D*b;
. P1 C' g' e  ~/ I7 o4 bJ=B*u0+f;6 @% D7 }, ]7 ?# Z" s
while norm(J-u0)>=eps% ~- X8 B$ b9 ^  ^* o
x0=J;
" e9 j! N. g: \, M6 Q% `6 m3 [J=B*u0+f;/ B. g2 p, t% N$ U
end
7 I2 ?, E( t; v8 ?: Ireturn
作者: 舒米牛牛    时间: 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