[热点讨论]请大家帮我看下我的MATLAB卷积分问题出在哪里？谢谢！ [复制链接]

用MATLAB自带卷积分函数计算弹性半空间在一三角脉冲荷载下自由表面的竖向位移，与文献正确结果相差10的5次方倍，但自己编写卷积分代码，反而得到文献结果。请大家帮我看下问题出在哪？谢谢！ k "'q   
qYGnebn@\  
另外，为何MATLAB自带卷积分函数结果矩阵的维数会是被卷积两矩阵维数之和减1？ 6./h0kD`  
X./8 PK?&  
1. 问题描述，如下图1所示： r-H~MisL  
z!D >l  
Uk6Y6mU V  
Z\6azhbI}  
2. 格林函数，如下图2所示： :sk7`7v  
(ZT*EFhb(  
%:YON,1b=7  
('OPW&fRG  
3. 参考文献正确解，如下图3所示： ~EpMO]I  
leg@ia  
vn6/H8  
V0c*M>V  
4. 本人直接用自编卷积分MATLAB代码算的出解（如下图4）及相应代码： "Yq-s$
yBi  
3)EslBA7i  
q~_Nv5r%O  
&CCB;Oi%  
clear; CZxQz  
clc; DUM,dFIlvF  
AL&}WbUC  
%loading history >.\G/'\?  
r/Qq-1E  
dt=0.01;
Ep.,2H  
ti=0.0008; {~[H"h537t  
te=12; 'nBJ[$2^  
t=ti:dt:te; %4i
mlP  
pt=(10*t./1.5).*(t>=0 & t<=1.5)+(20-10*t./1.5).*(t>1.5 & t<=3)+0.*(t>3); IP-
CN  
m=length(pt); /vD5C  
xJJlVP  
%load distribution in space #:DDx5%x<b  
UVxE~801Y  
rp=0.1; n|yl3v  
ri=0; [ C]=p  
rc=0.1; 1Jd82N\'  
dr=0.001; y%v<Cp@R  
r=ri:dr:rc; rV[#4,}PF  
pr=1*(r<=rp)+0*(r>rp); zPH1{|H+l  
n=length(pr); )PB&w%J  
uy~5!i&  
%load function with respect to t and r {KdC51"Nv  
Y%anR|  
p=pr.'*pt; l/DV ?27  
`m`jX|`  
%green's function s7D_fv4e  
T5Sa9\`>  
G=1; _+wv3? c"  
cs=1; [/6$P[  
for i=1:1:m R]m`v: 9  
  for j=1:1:n JrOp-ug  
    u(j,i)=heaviside(cs*t(i)-r(j))/(pi*G*sqrt(t(i)^2-(r(j)/cs)^2)); &
gIu<*u<  
  end f/eT4y  
end =}Bq"m  
Gxy>aS3  
%convolution and response of displacement 7.hVbjy'-  
D2[uex  
for i=1:1:m }8fxCW*|  
  for j=1:1:n )wCA8  
    v(j,i)=0; N@58R9P<p  
    for k=1:1:i B#B$w_z  
      for g=1:1:j tA3]6SIK@  
        v(j,i)=v(j,i)+p(g,k)*u(j-g+1,i-k+1)*dr*dt; aVO5zR./)  
      end; :_p3nb[r  
    end; #fzvK+  
  end; Haaungb"  
end; ;7K5Bo  
<@A/`3_O)  
%plot the response history QKE$>G  
vObP(@0AM  
plot(t,v(n,); nwOr  
xlabel('t/s'); w`c9_V  
ylabel('v/m'); |hiYV  
grid on; p!
zC  
title('response of point A or C'); Iwize,J~X  
D7 [n^WtL  
5. 直接用MATLAB自带卷积分算的出解（如下图5）及相应代码： 9K Ih}Q@P  
hG2btmBh
t  
gmj a2F,  
|\XjA4j  
clear; c zL[W2l  
clc; 4ME8NEE  
jf$6{zO6j  
%loading history  C!Y|k.`p  
JUF[Y^C  
dt=0.01; V4l`Alr\L  
ti=0.0008; ~ifq_Ag.  
te=12; [WRs1$5  
t=ti:dt:te; o-lb/=K+  
pt=(10*t./1.5).*(t>=0 & t<=1.5)+(20-10*t./1.5).*(t>1.5 & t<=3)+0.*(t>3); [j 'Og
m7"  
m=length(pt); 9c=Y+=<  
Zo-E0[9  
%load distribution in space 8}{';k  
}rn}r4_a  
rp=0.1; `$nMTx]Y  
ri=0; N+"Y@X yg  
rc=0.1; yodJGGAzk  
dr=0.001; #CW]70H`  
r=ri:dr:rc; +r2E5s  
pr=1*(r<=rp)+0*(r>rp); C ])Q#!D|  
n=length(pr); 67<zBw2  
e ! 6SJ7xC  
%load function with respect to t and r 4)]g=-3  
V/"XC3/n*  
p=pr.'*pt; ~SUA.YuF  
]BO{Q+?d2  
%green function 0u'4kF!P!  
]6pxd \Q  
G=1; 'R*xg2!i  
cs=1; =yz#L@\!  
for i=1:1:m nAoGG0$5  
  for j=1:1:n 0{BPT>'  
    u(j,i)=heaviside(cs*t(i)-r(j))/(pi*G*sqrt(t(i)^2-(r(j)/cs)^2)); -l`@pklQ  
  end 5w@
;B  
end 4CW/  
DcQ^V4_  
%convolution and response of displacement U#Wc!QN-t  
&Sa<&2W4S  
v=conv2(u,p); uQ vW@Tt  
\Y Cj/tG8  
%plot the response history /21d%T:}  
w>IkC+.?  
rr=r(n); })M$#%(  
[tpu,rpu]=meshgrid(ti:dt:2*te-dt,ri:dr:2*rc); Q2Yv8q_}Uq  
[X,Y,Z]=meshgrid(linspace(min(tpu(),max(tpu()),linspace(min(rpu(),max(rpu()),linspace(min(v(),max(v())); |n}W^}S5  
V=Y; A AH-Dj|&l  
h=contourslice(X,Y,Z,V,tpu,rpu,v,[0 0]+rr); !_?<-f(  
set(h,'edgecolor','k'); HvKueTQ  
contourslice(X,Y,Z,V,tpu,rpu,v,[0 0]+rr); o~NeS|a  
xlabel('t/s'); MPAZ%<gmD  
ylabel('r/m'); U^+9l?ol  
zlabel('v/m'); MN$j{+!Q  
axis([0 12 0 0.1 0 12e4]); f0IljY!.  
view(0,0); ^;6~=@#*C  
grid on; d?v#gW
 
title('displacement response history of point A or C');