EDA365电子论坛网
标题:
一维抛物有限元编程，检查半天也不知道哪错了，大家可以帮忙看看吗，谢谢了
[打印本页]
作者:
justlikethisis
时间:
2020-9-11 14:12
标题:
一维抛物有限元编程，检查半天也不知道哪错了，大家可以帮忙看看吗，谢谢了
function L2=lunwen(N,M)
( ^) Y( `! \) n, W- a7 h! C1 ?
lx=0;rx=1;lt=0;rt=1;
) d: V4 d* V* H& n$ u
h=(rx-lx)/N;
2 }0 U+ }/ w' K9 n! Y
th=(rt-lt)/M;
3 j% Z$ a* k* V* J( w- d
x=0:h:1;
9 A% b/ x% W4 B# N) S% @
t=0:th:1;
1 R. c- O4 D$ g$ I. j3 k
A=zeros(N+1,N+1);
7 X. }1 K4 K* M+ A$ h p4 K7 Y
B=zeros(N+1,N+1);
. F% F" l- A5 U& |& S* r
D=zeros(N+1,N+1);
* r; g9 G$ f! ], y% i& @
E=zeros(2*N+2,2*N+2);
* @" M- P1 F1 o2 O( ~2 T6 p
F=zeros(N+1,1);
: r7 F4 u: B+ s0 `) b0 _
u3=zeros(N+1,1);
% C! H1 ^2 o& ^" D& ?& o
P=zeros(2*N+2,1);
! ^, A* n, _8 g) m
X=zeros(2*N+2,1);
# a' `- q4 c: {: U2 L- Y
X1=zeros(2*N+2,1);
% x: m, ]# ]3 c* K% q
X2=zeros(2*N+2,M+1);
( A0 v# U _4 h4 X, z/ X7 Q
for i=0:N
4 ~' K5 K6 P, ^* t0 Z m6 B* L
u3(i+1,1)=sin(pi*i*h);
) p- ?4 K5 C' j' Q( q
end
2 Q1 }( k8 ~. P% o/ W, ~* F
X(N+2:2*N+2,1)=u3(1:N+1,1);
" T0 o6 c; d' f& x+ @8 `
for i=1:N
1 s- p @; j1 N5 T" M
e=[x(i),x(i+1)];
7 v* S' k/ E3 t( E2 E
A([i,i+1],[i,i+1])=A([i,i+1],[i,i+1])+ganddu1(e);
_: b. @2 {/ i3 @7 @
B([i,i+1],[i,i+1])=B([i,i+1],[i,i+1])+gangdu22(e);
9 N# b! p f; R
D([i,i+1],[i,i+1])=D([i,i+1],[i,i+1])+gangdu33(e);
" p/ p, x7 O$ b: I3 @7 ^: }
end
' D1 F4 R X, v( L4 G/ O' {$ k
% E(1:N+1,1:N+1)=A;
p0 l' ^3 R8 g- _; a1 u7 o7 ]
% E(1:N+1,N+2:2*N+2)=-B;
: G) B4 b% E5 d% q
% E(N+2:2*N+2,1:N+1)=-th*B;
9 ^& n4 w1 N7 }# Z
% E(N+2:2*N+2,N+2:2*N+2)=B-th*B+th*D;
: L! ]& h. y- |% E$ Y5 W5 ]5 t/ [* ~
E=[A,-B;-th.*B,B-th.*B+th.*D];
! K2 E4 i/ C2 D
for j=1:M
2 J4 j! }4 h* \, a0 n2 c
F=zeros(N+1,1); %%记得归零
6 o7 z" A; ]4 [% h+ M' I& j; C
y=exp(-t(j+1));
; a5 a; m8 i* [! c8 }. i9 o
for i=1:N
9 i6 |: W/ B+ X5 S2 O, L" n7 S' Y
e=[x(i),x(i+1)];
4 ~% k5 T2 T; N- g+ c6 e. i
F([i,i+1],1)=F([i,i+1],1)+y*Hezai3(e);
! r4 f9 J' \% i' d$ _
end
) q! h# J+ O# z$ G; ~5 W. v2 w ]
X(1:2.*N+2,j+1)=E\[zeros(N+1,1);-th.*F+B.*X(N+2:2.*N+2,j)];
3 A" H/ W( I! E3 x3 o
end
7 ~, s! n" x/ W4 V8 k
X1=[A\B*X(N+2:2*N+2,1);X(N+2:2*N+2,1)];
' g4 V. s, Y: o" |1 i4 i1 I8 N _
X2=[X1;X];
2 f( S6 ]3 s% R% v
错误使用 vertcat
. w5 b9 d) n% q7 J1 o
要串联的数组的维度不一致。
: |4 m( U1 H$ K
$ Q' Y: P! z' k7 g7 |; d
出错 lunwen (line 39)
* ]! g5 G* `& C5 t% r) Y: F
X(1:2.*N+2,j+1)=E\[zeros(N+1,1);-th.*F+B.*X(N+2:2.*N+2,j)];
" g6 Q( G' e# j
检查半天也不知道哪错了
Y4 R$ ?) f) I. Z* P' h/ ?, V
9 k4 E+ ]2 V& f6 O4 J
, _1 P3 A+ d/ Z3 b! s7 N: t( y
作者:
qq666888qqw
时间:
2020-9-11 15:11
你的程序没传完吧
作者:
somethingabc
时间:
2020-9-11 16:35
程序好像没传完
欢迎光临 EDA365电子论坛网 (https://bbs.eda365.com/)
Powered by Discuz! X3.2