数学建模社区-数学中国

标题: 求助：关于分块矩阵的还原，为什么实行不了 [打印本页]
作者: 舒米牛牛    时间: 2009-12-17 20:17
标题: 求助：关于分块矩阵的还原，为什么实行不了
%运用Jacobi迭代拟合出来的u关于x和y的矩阵
: p; c! W( S7 Z5 A# p& bfunction J=jac(A,b,u0,eps)- B! u3 K0 A/ ~
if nargin==37 I3 ]! {$ h0 ]  N
    eps=1.0e-8
  Q- {- q7 X  q/ n4 kelseif nargin<3
" |% m# O& [/ s# @    'error'- G; ]# n! l+ K) F6 _% p( T3 V3 m
    return
7 `% Q: ]8 w  j  ]1 x+ _& Gend
. G0 h. s6 k+ c/ V# O. T$ w0 {6 {- S  M4 Q
%定义内部节点矩阵u0/ @$ R* M0 g9 @
h=1;k=1;* r4 A0 w7 i: ^
x=0:h:17;y=0:k:10;
5 \1 r, {# ]6 R( ke=length(x)-2;f=length(y)-2;
$ b! w" ^! f9 @3 z2 vu=zeros(e,f);
7 z# M9 N7 T5 H: V/ R. Nu0=u;" x9 J3 t6 l2 J' ]* r$ H1 v

$ I5 B2 V- D5 r7 x5 ~%定义外部节点p. L2 e3 U7 z) |+ g& G$ Q
p=zeros(e+2,e+2);6 p3 u- g" `! l8 L& F: Z( H2 a
p(1,1:f+2)=0;p(e+2,1:f+2)=0;       ; W2 x# t, X( ~& y! H( K
p(1:e+2,1)=100;p(1:e+2,f+2)=100;& u8 R. m. Y) o/ u  q5 g9 z

% J- V! h# A2 k5 I) `3 T%定义系数矩阵A. B7 J% r" X6 n; g0 r. B1 U* s9 N
A=zeros(e*f,e*f);
/ X* P' T. \+ z+ e! n7 [* z3 XB=mat2cell(A,ones(e*f/e,1)*e,ones(e*f/e,1)*e);
4 g( t. Z, q( w+ i2 td1=ones(e,1);d2=ones(e-1,1);" R: ]0 ]0 R+ x) y/ k# s# o0 W
M=4*diag(d1)-diag(d2,1)-diag(d2,-1);
. E* S  H0 {+ c: pN=-eye(e);
' P# _  Q  }. j0 Z5 h3 b' ^B{1,1}=M;B{e,e}=M( j) O5 ^. a. g4 y3 S; N& Z
for i=2:e-1
5 a1 D1 ?0 _! k4 w, ?    B{i,i}=M;B{i-1,i}=N;B{i+1,i}=N" V( z. [6 E: L6 p) K- X; o
end7 Q- P+ S/ h$ [/ e4 f5 I
A=cell2mat(B);  & e1 T, _/ e; ~1 K8 ^0 H+ b
这里总是显示% o5 |+ G$ D* l
??? function J=jac(A,b,u0,eps)0 q, U, Q- M; c) q1 w- {
    |
1 _  J: L9 }+ s+ J5 m# @( d: \Error: Function definitions are not permitted at the prompt or in scripts., d4 }2 \2 r" Y! [& @8 ^0 O
3 H. r% C  J8 Q& i* N
%定义b+ D6 O0 a% Y4 y
b=ones(e*f,1);
: d' b' N7 q; t4 T4 Efor i=1:e
& ]+ Z3 _+ E- h# o, O5 k7 |    for j=1:f
3 ~7 w- ^% q$ ~* B2 {: w2 {        b(i+j)=p(i,j+1)+p(i+1,j)+p(i+2,j+1)+p(i+1,j+2)  I8 t8 p+ t" T+ ^- {1 c
    end% {  v% q* r% y, X1 s" w3 D
end0 _# E$ Z% H# B. K- Y3 f$ D4 @, J
%运用Jacobi迭代法计算
. H6 M5 U9 M: v4 o# `5 fD=diag(diag(A));
( {  `, m( n* c2 i4 N, u/ o+ P* BD=inv(D);( `5 I+ W9 G: F0 N4 g
L=tril(A,-1);
* |! }* d/ y* `) V/ p+ tU=triu(A,1);
+ x. r6 j$ t/ a$ O1 cB=-D*(L+U);4 g* w0 F# u/ A8 b+ X9 L( O* ~
f=D*b;
& F8 B- J7 x9 O  W6 OJ=B*u0+f;7 }1 {( j7 G% S
while norm(J-u0)>=eps
6 V% X" x4 O' lx0=J;
5 Y# W( X  Y/ }1 U- y! rJ=B*u0+f;
% v8 K3 P! j* send
" R$ x  L8 n$ N3 j& sreturn
作者: 舒米牛牛    时间: 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