EDA365电子论坛网
标题:
在赋值 A(:) = B 中，A 和 B 中的元素数目必须相同。
[打印本页]
作者:
purpose_857
时间:
2023-3-1 10:13
标题:
在赋值 A(:) = B 中，A 和 B 中的元素数目必须相同。
用matlab求三元一次微分，出现了如下问题，请大神帮忙看下：
/ v& C" h5 a. y2 Q7 k2 p
问题：在赋值 A(: ) = B 中，A 和 B 中的元素数目必须相同。
2 [! _! J/ x' g
& K. B/ e! P- @# }2 W& k' d
出错 longgekuta_lunwen2 (line 41)
8 V; b) p9 M1 q4 ~* R+ U; ^
y(i+1) =y(i)+h/6*(K1+2*K2+2*K3+K4);
, T8 _( A2 C0 m7 Z: Z
& ^- j; q& i1 G7 w9 ]$ K+ U
代码如下：clc;
0 u# }& X; N) b
clear all;
U4 ]% ]7 l- r1 V* v
global r midu E Es u c q t w i n m1 m2 A B U1 S1 o j %o为弹塑性界面速度（C)与缓冲层弹性波速（Cd）的比值 f为应变值
4 H/ u: Y6 [- e0 H; p$ Q( z" e
r=0.2; %落石半径
4 u. D% J( b0 k5 c- D/ u; X# X
midu=1640;%密度
5 @% ~# Q- I1 k3 I: L7 L
E=35e6; %弹性模量
+ Y. F" o2 @4 M" r+ p- q
Es=8e6; %压缩模量
5 l; F, g i& p) d/ P! Q
u=0.32; %泊松比
0 [. |% w7 I4 r: o% X* x& O
c=12e3; %粘聚力
* l- T* d! g+ s6 J# | v
q=30 ; %内摩擦角
6 ~9 e7 w7 U$ o( A" L U" l
t=sind(q);
% `; z; G( {6 ?$ g( n2 Z2 q" Z% h
w=3*sind(q)/sqrt(3+sind(q)^2); %屈服强度
- c/ P2 E: e* y
i=3*c*cosd(q)/sqrt(3+sind(q)^2);%剪切强度
) d6 Q& Z! W* j. h
& w; ]; {9 L, Y4 [ S. L: k9 e1 O
n=3/(3+2*w);
9 }) {+ J1 [/ @; `
m1=6*w/(3+2*w);
. y0 m/ A/ M( L- g4 |0 l
m2=6*i/(3+2*w*Es);
L7 k! p" m) Z
u/ h+ `: n5 s$ y% Q1 K! W+ D6 J& b# c
A=6*(1+u)*(1-2*u)*o^2*i/((3+2*w)*(1-o)*E*((1-2*u)*(3-6*w/(3+2*w))*(1+o)-(6*w*o^2*(1+u))/(3+2*w)));
& b0 |2 s2 g( r
B=-3/2*A;
( S( G9 ~' g8 E! |6 z- y# ?9 {1 u
U1=3*(1-o^2)*A;
7 O7 c& ` L9 z) y3 f+ x9 U7 U7 J
S1=-E/(Es*(1+u)*(1-2*u))*((1+u)*A*o+2*(1-2*u)*B/(3*o^2)+2*u*B);
4 s& ^+ Z1 c* s' v; x2 |
. A# w( y2 q$ {! q
h=j-1;
" u2 e# ?" k/ f: `
N =10;
2 M5 l% u" Q0 a4 l1 C* w
x1=1; %初始值1
& M6 {; c0 ` m
y1=U1; %初始值1
' S6 _1 n* X4 b P e9 V5 e
z1=S1; %初始值1
4 O1 {0 F) m$ T, Y# G( E0 W$ z2 k
x2=j; %初始值2
: ~( E+ y) f% X! t7 V" |1 e, \
y2=j; %初始值2
3 g. C$ J4 f: w, [! E9 g) G5 n
x = ones(1,N);
- x. w9 D+ ~0 A
y = ones(1,N)*0;
# q* @4 X) o3 j$ L& a" X4 {& P5 o
z = ones(1,N)*0;
8 D+ u9 L& X5 t" N2 [6 x
5 u- b0 b- G V& L
for i=1:N
$ [9 g7 s; P! D
[K1,L1] = F(1,y(i),z(i))
" F3 [/ f) |# W
[K2,L2] = F((j+1)/2,y(i)+h*K1/2,z(i)+h*L1/2)
- _- e$ m0 Y+ U/ P
[K3,L3] = F((j+1)/2,y(i)+h*K2/2,z(i)+h*L2/2)
* J8 |/ B% A2 {* s7 J* O
[K4,L4] = F(j,y(i)+h*K3,z(i)+h*L3)
; M; t3 z1 f" M0 @: w$ u
4 O4 h6 Z5 G5 X/ U+ q0 C
y(i+1) =y(i)+h/6*(K1+2*K2+2*K3+K4);
3 `8 i7 s8 A. g0 c" [# d* I; Q
z(i+1)= z(i)+h/6*(L1+2*L2+2*L3+L4);
6 N6 n l4 W. Y8 P
end
2 l; a# }8 E# N& E% F
1 f' T. Y. C. \- J4 X
figure(1);
- n) _3 Q) C0 _ a/ v( w
comet3(x,y,z,0.1)
$ r) B; h" k% M" B+ D4 {
function [K,L] = F(X,U,S)
6 V" c! n w6 R& W
global r midu E Es u c q t w i n m1 m2 A B U1 S1 C %C为弹塑性界面速度
( d6 S" q+ ~: u, r2 ^" ?. _
K = (2*U/X+3*(m1*S+m2)*(X-U)/(X*(n-3*S+2*i/Es)))./((3*n*midu*C^2)*(X-U)^2/(Es*(n-3*S+2*i/Es)^2)-1);
: \/ r: y2 {& v U
L = ((2*U/X+3*(m1*S+m2)*(X-U)/(X*(n-3*S+2*i/Es)))./((3*n*midu*C^2)*(X-U)^2/(Es*(n-3*S+2*i/Es)^2)-1))*n*midu*C^2*(X-U)/(Es*(n-3*S+2*i/Es))-(m1*S+m2)/X;
( \2 N+ i1 Y/ {3 J
end
1 S' q8 N; R. w; E! c* ~, _
请大神赐教，感谢!
% z& C+ U* g; k3 t5 h
作者:
whatever_
时间:
2023-3-1 11:27
仅供参考，主程序里的j，o，C三个变量未定义，且最后一个变量未在主程序里声明为全局变量
作者:
purpose_857
时间:
2023-3-1 11:38
whatever_ 发表于 2023-3-1 11:27
- q) ^% A9 c% j4 k3 R" c$ ?
仅供参考，主程序里的j，o，C三个变量未定义，且最后一个变量未在主程序里声明为全局变量
- w$ X3 d7 i1 B7 l2 F% C0 V
这里 o, C 是未知量，给 j 赋初始值的情况下，这里微分能求出来吗？
& [' t0 j6 Y! r2 w; a; L3 {
作者:
peerless2021
时间:
2023-3-1 13:49
只是单纯从代码角度出发，你这里的o, C需要赋值，否则就会导致空集，仅供参考
欢迎光临 EDA365电子论坛网 (https://bbs.eda365.com/)
Powered by Discuz! X3.2