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

用MATLAB自带卷积分函数计算弹性半空间在一三角脉冲荷载下自由表面的竖向位移，与文献正确结果相差10的5次方倍，但自己编写卷积分代码，反而得到文献结果。请大家帮我看下问题出在哪？谢谢！ Z/7dg-$?'0  
I="oxf#q  
另外，为何MATLAB自带卷积分函数结果矩阵的维数会是被卷积两矩阵维数之和减1？ 0;<OYbm3<  
a_{6Qdl  
1. 问题描述，如下图1所示： {*$9,  
a:b^!H>#  
?:/|d\,7@  
GS4_j
vD-  
2. 格林函数，如下图2所示： jA<T p}$!  
Egf^H>,.M  
\8>oJR 6  
{R8=}Qo  
3. 参考文献正确解，如下图3所示： 6c&Y  
fGTOIi@#  
Yf=FeH7"  
HY*\ k#  
4. 本人直接用自编卷积分MATLAB代码算的出解（如下图4）及相应代码： UJn/s;$.e  
-TS? fne)  
8gI\zgS  
nvH|Ngg Q  
clear; hfv%,,e  
clc; JiA'BEJN  
/WYh[XKe  
%loading history v)+@XU2wZ  
 D%gGRA  
dt=0.01; H(&Z:{L  
ti=0.0008; 'F7VM?HBfg  
te=12; 8(Fu  
t=ti:dt:te; 11{y}J  
pt=(10*t./1.5).*(t>=0 & t<=1.5)+(20-10*t./1.5).*(t>1.5 & t<=3)+0.*(t>3); NnOI:X {  
m=length(pt); vYdlSe=6G  
Dft%ip2  
%load distribution in space L {qJ-ln:  
lkwh'@s.  
rp=0.1; H;y}-=J+  
ri=0; {g_@Tuu  
rc=0.1; *Ru2:}?MpS  
dr=0.001; hDvpOIUL1  
r=ri:dr:rc; %E.S[cf%8&  
pr=1*(r<=rp)+0*(r>rp); Gkmsaf>  
n=length(pr); >|nt2  
"lrA%~3%[P  
%load function with respect to t and r V.2[ F|P;3  
l;0y-m1  
p=pr.'*pt; ]7vf#1i<  
_Ex|f5+  
%green's function 7=3O^=Q^Q  
7xT[<?,  
G=1;
%Rarr  
cs=1; Ow)R|/e/  
for i=1:1:m l"5y?jT  
  for j=1:1:n nh0&'hA  
    u(j,i)=heaviside(cs*t(i)-r(j))/(pi*G*sqrt(t(i)^2-(r(j)/cs)^2)); )5GQJiY  
  end agT7=hX].  
end 1.0J2nZpt  
Q7(eq0na  
%convolution and response of displacement {i;6vRr  
CjKRP;5  
for i=1:1:m Y&GuDLUF  
  for j=1:1:n TGpSulg7  
    v(j,i)=0; ,C:o`fQ\  
    for k=1:1:i W_}/O'l{  
      for g=1:1:j Y
1y E  
        v(j,i)=v(j,i)+p(g,k)*u(j-g+1,i-k+1)*dr*dt; 4U{m7[  
      end; .CS v|:'1  
    end; +*.1}r&  
  end; g`3H(PVg  
end; &O*ENpF  
&h(g$-l?[  
%plot the response history ]! )xr  
DY.58IHg1  
plot(t,v(n,); "i%jQL'.  
xlabel('t/s'); l{Er+)a  
ylabel('v/m'); C0(s
AF@  
grid on; sUc
iFAb  
title('response of point A or C'); ET+'Pj3  
'hI
U_  
5. 直接用MATLAB自带卷积分算的出解（如下图5）及相应代码： iaRR5
D-  
kFwxK"n@C  
%w:'!X><  
9|3o<  
clear; t3>$|}O]t  
clc; Z Xb}R^O-  
=:/>6H1x  
%loading history P^zy;Qs7  
L$hc,  
dt=0.01; A{(T'/~"  
ti=0.0008; q~h:<,5  
te=12; >qpqQ; bm  
t=ti:dt:te; Mpm#GdT  
pt=(10*t./1.5).*(t>=0 & t<=1.5)+(20-10*t./1.5).*(t>1.5 & t<=3)+0.*(t>3); 8Zw]f-5x\  
m=length(pt); \O? u*  
;($1Z7j+  
%load distribution in space |_nC6;  
wT/6aJoX  
rp=0.1; +nQ!4  
ri=0; u>o<tw%Y  
rc=0.1; }p{;^B  
dr=0.001; zt?H~0$LB  
r=ri:dr:rc; +'%\Pr(  
pr=1*(r<=rp)+0*(r>rp); ^1VbH3M  
n=length(pr); 1Is%]6  
e1uMR-Q  
%load function with respect to t and r GA@ Ue9  
sOQcx\dK  
p=pr.'*pt; nq@5j0fK
 
M=[th  
%green function 5#!ogKQ(i  
!g2a|g  
G=1; o(Kcs-W2  
cs=1; oyW00]ka  
for i=1:1:m 9-93aC.|}  
  for j=1:1:n &^+3errO  
    u(j,i)=heaviside(cs*t(i)-r(j))/(pi*G*sqrt(t(i)^2-(r(j)/cs)^2)); Abf1"#YImy  
  end Spo+@G  
end uP6-cs  
L|J~9FM  
%convolution and response of displacement TPK@*9rI  
+*D4(  
v=conv2(u,p); ?gG,
t4D  

F[]&1  
%plot the response history MD4\QNUa)*  
S&01SX6  
rr=r(n); ^@"c`  
[tpu,rpu]=meshgrid(ti:dt:2*te-dt,ri:dr:2*rc); `Cg^in\  
[X,Y,Z]=meshgrid(linspace(min(tpu(),max(tpu()),linspace(min(rpu(),max(rpu()),linspace(min(v(),max(v())); "
K ?#,
_  
V=Y; 85hQk+Bu4  
h=contourslice(X,Y,Z,V,tpu,rpu,v,[0 0]+rr); n$W"=Z;`  
set(h,'edgecolor','k'); 0x71%=4H^x  
contourslice(X,Y,Z,V,tpu,rpu,v,[0 0]+rr); U`1l8'W}:#  
xlabel('t/s'); &CUC{t$VHX  
ylabel('r/m'); 74]a/'
4  
zlabel('v/m'); JY@X2'>v/  
axis([0 12 0 0.1 0 12e4]); @d)LRw.I  
view(0,0); g@u;Y5  
grid on; N&x:K+Zm.  
title('displacement response history of point A or C');