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

用MATLAB自带卷积分函数计算弹性半空间在一三角脉冲荷载下自由表面的竖向位移，与文献正确结果相差10的5次方倍，但自己编写卷积分代码，反而得到文献结果。请大家帮我看下问题出在哪？谢谢！ Og`6>?>97  
zL@ZNH  
另外，为何MATLAB自带卷积分函数结果矩阵的维数会是被卷积两矩阵维数之和减1？ dpO ZqhRs.  
RNyw`>  
1. 问题描述，如下图1所示： T(<C8  
zkdyfl5  
)w8h2=l  
iBy:HH  
2. 格林函数，如下图2所示： ,H3~mq]  
4Lx#5}P  
6'e}!O  
zY<=r.m4  
3. 参考文献正确解，如下图3所示： 0R,?$qM\  
c}II
"P  

VP$`.y  
fjUyx:  
4. 本人直接用自编卷积分MATLAB代码算的出解（如下图4）及相应代码： +h6cAqm]  
^/wvHu[#  
05zBB  
9HlRf6S  
clear; R4<}kA,.  
clc; ,y[wS5li  
F6gboo)SD  
%loading history +8FlDiP  
9L}=xX`>?  
dt=0.01; '3f"#fF6  
ti=0.0008; i#t)tM
"  
te=12; ]@W.5!5H  
t=ti:dt:te; ,%+i}H,3  
pt=(10*t./1.5).*(t>=0 & t<=1.5)+(20-10*t./1.5).*(t>1.5 & t<=3)+0.*(t>3); hsVf/%  
m=length(pt); o;XzJ#P  
g/b_\__A  
%load distribution in space JDi|]JY  
-d+q+l>0  
rp=0.1; 9PA\Eo|Yb  
ri=0; Qwn/ ,  
rc=0.1; t)n!];  
dr=0.001; qV$\.T>x  
r=ri:dr:rc; eI@LVi6<b  
pr=1*(r<=rp)+0*(r>rp); fA u^%jiU  
n=length(pr); Y~T;{&wi  
\: H&.VQ"  
%load function with respect to t and r K.cMuh  
"CdL?(  
p=pr.'*pt; \N\Jny  
ic:_v?k  
%green's function DiyviH  
VRYj&s'@  
G=1; -hf)%o$  
cs=1; .17WF\1HC.  
for i=1:1:m !"2nL%PW~  
  for j=1:1:n -{i;!XE$SR  
    u(j,i)=heaviside(cs*t(i)-r(j))/(pi*G*sqrt(t(i)^2-(r(j)/cs)^2)); "-Lbz)k  
  end 5MH\Gqe7  
end ZM/*cA!"  
^+
zF;Q'  
%convolution and response of displacement n|vIo)  
}fW@8ji\  
for i=1:1:m -X~VXeg  
  for j=1:1:n P1b5=/}:V  
    v(j,i)=0; <im BFw  
    for k=1:1:i 1&~u:RUXe  
      for g=1:1:j vdloh ,  
        v(j,i)=v(j,i)+p(g,k)*u(j-g+1,i-k+1)*dr*dt; /jQW4eW0  
      end; [q/=%8qLUA  
    end; /O1r=lv3Z  
  end; dnVl;L8L3  
end; (yv&&Jc  
!5Ko^:+Y  
%plot the response history @Y'BqDFlZ  
W8Z&J18AU  
plot(t,v(n,); DUc -D==  
xlabel('t/s'); $UMxO`F  
ylabel('v/m'); CPVR  
grid on; ~3Y)o|D3  
title('response of point A or C'); 9 veq  
UdmYS3zs  
5. 直接用MATLAB自带卷积分算的出解（如下图5）及相应代码： 7hq*+e  
gG0P &9xz  
4tTK5`7N  
#`5{?2gS9  
clear; LB_ylfg  
clc; lzz rzx^  
k&4@$;Ap  
%loading history +(iM]L$Fw%  
>|L,9lR_b  
dt=0.01; 12*'rU;*  
ti=0.0008; oHkF>B [  
te=12; TY[1jW~{r  
t=ti:dt:te; 
U+t|wK  
pt=(10*t./1.5).*(t>=0 & t<=1.5)+(20-10*t./1.5).*(t>1.5 & t<=3)+0.*(t>3); d]0.6T1[K  
m=length(pt); Gxu&o%x[  
q;a`*gX^  
%load distribution in space R9o3T)9V  
*hVW>{a  
rp=0.1; P SDzs\s  
ri=0; k:nR'TI  
rc=0.1; " d3pkY  
dr=0.001; ;7"}I  
r=ri:dr:rc; |:SBkM,  
pr=1*(r<=rp)+0*(r>rp); 7@u0;5p|  
n=length(pr); klduJT >  
Yxye?R-:  
%load function with respect to t and r {O,D9
<  
V-U,3=C  
p=pr.'*pt; pOlo_na}[  
>OVi{NyT  
%green function )A9K
9pZj  
$v?+X20  
G=1; )cnB>Qul  
cs=1; 9WsPBzi"T  
for i=1:1:m 5|!x0H;  
  for j=1:1:n $d M: 5y  
    u(j,i)=heaviside(cs*t(i)-r(j))/(pi*G*sqrt(t(i)^2-(r(j)/cs)^2)); TTaSg\K  
  end s5T$>+ a  
end 6LRI~*F=3  
`#&pB0.y  
%convolution and response of displacement tuuc9H4B  
&B\tcF  
v=conv2(u,p); E.]sX_X?  
F gM<2$h  
%plot the response history 7pDov@K<{  
Q}~of}h/  
rr=r(n); je1f\N45  
[tpu,rpu]=meshgrid(ti:dt:2*te-dt,ri:dr:2*rc); LuQ4TT  
[X,Y,Z]=meshgrid(linspace(min(tpu(),max(tpu()),linspace(min(rpu(),max(rpu()),linspace(min(v(),max(v())); 0x!XE|7I  
V=Y; & IVwm"  
h=contourslice(X,Y,Z,V,tpu,rpu,v,[0 0]+rr); /238pg~Cw5  
set(h,'edgecolor','k'); m5lMh14E  
contourslice(X,Y,Z,V,tpu,rpu,v,[0 0]+rr); c#M'Mye  
xlabel('t/s'); sw{,l"]<  
ylabel('r/m'); [L 0`B9TD~  
zlabel('v/m'); 8;YeEW5  
axis([0 12 0 0.1 0 12e4]); [6Y6{
.%~  
view(0,0); 6'UtB!gr  
grid on; Be+CV">2  
title('displacement response history of point A or C');