数学建模社区-数学中国
标题:
求助:关于分块矩阵的还原,为什么实行不了
[打印本页]
作者:
舒米牛牛
时间:
2009-12-17 20:17
标题:
求助:关于分块矩阵的还原,为什么实行不了
%运用Jacobi迭代拟合出来的u关于x和y的矩阵
0 G6 `/ e @' Q; D3 T" i
function J=jac(A,b,u0,eps)
0 U1 K: v+ J7 E9 H( L
if nargin==3
1 b6 d: z# ?! r( G5 [* r8 n1 H
eps=1.0e-8
' A: f; j7 P' |+ K- M \: V- G
elseif 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
%定义内部节点矩阵u0
9 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# Y
e=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 A
7 Y6 T" y/ r4 r- a* k% b
%定义外部节点p
6 \/ y; M$ _- U% n* c/ X+ l* |) ^
p=zeros(e+2,e+2);
! |) z; k" U# k5 m, T
p(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. Y
M=4*diag(d1)-diag(d2,1)-diag(d2,-1);
. e7 `) O: d& Q Y( c
N=-eye(e);
$ K7 [2 I: D% [* J
B{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- {- q
end
( 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 L
Error: 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 T
for i=1:e
5 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
end
1 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 w
U=triu(A,1);
0 h- Y, Q# F9 c q0 q7 h) D7 F8 u
B=-D*(L+U);
7 E* Y/ n: ]$ [1 q
f=D*b;
7 {1 i) N5 |" p3 D
J=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# B
x0=J;
- _% L( B5 v2 A, C3 P) y2 N" A! U
J=B*u0+f;
/ T. G- R! T2 }4 u+ u1 R; c# A. w" X
end
* 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