QQ登录

只需要一步,快速开始

 注册地址  找回密码
查看: 3042|回复: 3
打印 上一主题 下一主题

求助:关于分块矩阵的还原,为什么实行不了

[复制链接]
字体大小: 正常 放大

1

主题

2

听众

35

积分

升级  31.58%

该用户从未签到

自我介绍
20100103重要的日子
跳转到指定楼层
1#
发表于 2009-12-17 20:17 |只看该作者 |倒序浏览
|招呼Ta 关注Ta
%运用Jacobi迭代拟合出来的u关于x和y的矩阵
. A+ T" W0 M4 U7 zfunction J=jac(A,b,u0,eps)
* U& T. J* ~8 J0 i' }if nargin==3$ O6 c) u) p; k! ~
    eps=1.0e-8
1 `2 j* X$ O5 d9 a( W& n# e* \elseif nargin<3
  I' _1 \) y6 R, E: p7 u    'error'
7 X+ ^0 K6 L+ ^7 h/ o0 n    return
, W6 w, h7 x/ t- Xend
0 f5 u- ?0 S9 H" R6 D' Y( Q: D% p# ?1 R. ^# S
%定义内部节点矩阵u0
7 j. k; d# r4 v) Q: W* Y4 Uh=1;k=1;
! E9 h6 @1 L- ?. nx=0:h:17;y=0:k:10;4 n6 b9 ?7 S1 w* o
e=length(x)-2;f=length(y)-2;. S0 T1 O0 D! X: t: K0 S
u=zeros(e,f);
) \+ z2 Y! F! ~9 b% e( S0 r8 E% Du0=u;; _: y) f! D: _) T
: w( w/ X; g9 b# C* m# e0 L
%定义外部节点p1 v  n/ B1 I6 z7 e6 q+ [7 p
p=zeros(e+2,e+2);
, W4 D4 L, R: [+ Yp(1,1:f+2)=0;p(e+2,1:f+2)=0;       / }0 _- B7 b9 ?5 p
p(1:e+2,1)=100;p(1:e+2,f+2)=100;8 P6 g: b4 a+ T3 M
; v7 X( x9 q3 f7 u& m
%定义系数矩阵A) `3 r  C! I0 }3 [6 Z, U
A=zeros(e*f,e*f);
! P$ f- [/ o5 j/ ^B=mat2cell(A,ones(e*f/e,1)*e,ones(e*f/e,1)*e);3 s; l9 i8 j+ `: Z* J
d1=ones(e,1);d2=ones(e-1,1);
. _) Y( x& B" f6 y) }* ^. rM=4*diag(d1)-diag(d2,1)-diag(d2,-1);
6 W& u" e! y' E5 Q, f' S8 W9 n% XN=-eye(e);
( Z) T) L3 }3 q9 sB{1,1}=M;B{e,e}=M
6 Z2 o% y/ Z3 g9 d$ K$ x& ^for i=2:e-1
3 B8 t% h8 H  o3 E) P4 m    B{i,i}=M;B{i-1,i}=N;B{i+1,i}=N
3 p8 c' ]3 ]! k: ]  R. dend
) {* f! ~2 a: IA=cell2mat(B);  
6 N7 }9 `4 S6 r; }5 v这里总是显示2 X* J0 B( ~5 a) e
??? function J=jac(A,b,u0,eps)
% z! E  E$ n" |) `! X    |- f! i; }- }( X+ O1 r
Error: Function definitions are not permitted at the prompt or in scripts.
! g, N8 S" P; X* Y, \! |! p2 ?
8 I  Y7 Z8 `0 U8 \0 R% x
%定义b
+ k+ I: i( t$ w0 M, p+ F3 j8 A" Fb=ones(e*f,1);
! _% P2 u0 ~$ N0 a  sfor i=1:e/ E' n* W% E7 P% o9 B+ P8 Y
    for j=1:f
: \# m8 C6 V  A        b(i+j)=p(i,j+1)+p(i+1,j)+p(i+2,j+1)+p(i+1,j+2)
' Q3 |7 w3 d9 p1 \: o* E& s3 A: j    end% C1 D: l  h, }# @$ |, ?* A
end
* B: X0 i* G# ?2 O7 N) O%运用Jacobi迭代法计算) b2 y" [6 n+ X% h1 k
D=diag(diag(A));
$ S5 i. w1 ?6 y: OD=inv(D);
' j5 k4 i: g( QL=tril(A,-1);
7 w; F. x0 m* [, D1 |U=triu(A,1);8 t( u0 X& `$ k; c5 W2 b$ [
B=-D*(L+U);9 L) h& g/ s: h3 `
f=D*b;
0 F1 W% K; q3 X% lJ=B*u0+f;
) M( [1 a1 u# w' k9 k' X% Mwhile norm(J-u0)>=eps
7 w, o# r: [8 u8 l* O! F  {7 j0 fx0=J;  c+ ]% s) f( T
J=B*u0+f;  `' R% u+ Z2 o2 ~( J; H8 ?
end
/ t) h, g# f! o$ N0 }% D, c, Breturn
zan
转播转播0 分享淘帖0 分享分享0 收藏收藏0 支持支持0 反对反对0 微信微信

1

主题

2

听众

35

积分

升级  31.58%

该用户从未签到

自我介绍
20100103重要的日子
自己顶一下,拜托哪位高手指点一下,纠结这个矩阵的还原,想了好多方法还是不行~~实在想不出哪里出错了
回复

使用道具 举报

BenCam 实名认证       

9

主题

6

听众

89

积分

该用户从未签到

自我介绍
200 字节以内
不支持自定义 Discuz! 代码
回复

使用道具 举报

madio        

3万

主题

1312

听众

5万

积分

  • TA的每日心情
    奋斗
    2024-7-1 22:21
  • 签到天数: 2014 天

    [LV.Master]伴坛终老

    自我介绍
    数学中国站长

    社区QQ达人 邮箱绑定达人 优秀斑竹奖 发帖功臣 风雨历程奖 新人进步奖 最具活力勋章

    群组数学建模培训课堂1

    群组数学中国美赛辅助报名

    群组Matlab讨论组

    群组2013认证赛A题讨论群组

    群组2013认证赛C题讨论群组

    回复

    使用道具 举报

    您需要登录后才可以回帖 登录 | 注册地址

    qq
    收缩
    • 电话咨询

    • 04714969085
    fastpost

    关于我们| 联系我们| 诚征英才| 对外合作| 产品服务| QQ

    手机版|Archiver| |繁體中文 手机客户端  

    蒙公网安备 15010502000194号

    Powered by Discuz! X2.5   © 2001-2013 数学建模网-数学中国 ( 蒙ICP备14002410号-3 蒙BBS备-0002号 )     论坛法律顾问:王兆丰

    GMT+8, 2026-7-11 12:43 , Processed in 1.914761 second(s), 72 queries .

    回顶部