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

用MATLAB自带卷积分函数计算弹性半空间在一三角脉冲荷载下自由表面的竖向位移，与文献正确结果相差10的5次方倍，但自己编写卷积分代码，反而得到文献结果。请大家帮我看下问题出在哪？谢谢！ /(.mp<s0  
|bO"_U  
另外，为何MATLAB自带卷积分函数结果矩阵的维数会是被卷积两矩阵维数之和减1？ f)^_|8  
5 4L\Jx  
1. 问题描述，如下图1所示： Ljp%CI[i  
% a@>_  
w%JTTru  
i?
 K|TC`  
2. 格林函数，如下图2所示： }x07^4$j  
"HDcmIXg&  
@tZ&2RY1  
^h"`}[+  
3. 参考文献正确解，如下图3所示： ?'KL11@R  
@NNq z  
4UW_Do  
#0y)U;dA+w  
4. 本人直接用自编卷积分MATLAB代码算的出解（如下图4）及相应代码： XYZ4TeW\1  
+O*/"]h  
U'<KC"f:'!  
/Sc l#4bW  
clear; t+1 %RyKFB  
clc; TjwBv6h  
&C"L  
%loading history Y]B)'[=h  
JNJ=e,O,  
dt=0.01;  y"H*%]  
ti=0.0008; /Z@tv.f  
te=12; t3&LO~Ye  
t=ti:dt:te; $**r(HV  
pt=(10*t./1.5).*(t>=0 & t<=1.5)+(20-10*t./1.5).*(t>1.5 & t<=3)+0.*(t>3); Ljx(\Cm  
m=length(pt); 1ke g
9]  
&3TEfvz  
%load distribution in space ,I%g|'2  
+i@y@<l:+  
rp=0.1; <cqbUL  
ri=0; A*}.EClH  
rc=0.1; Dk(1}%0U/  
dr=0.001; >JC
  
r=ri:dr:rc; iC98_o_9  
pr=1*(r<=rp)+0*(r>rp); f;xkT  
n=length(pr); wv\w;'  
L;'"A#Pa  
%load function with respect to t and r ]y1OFKYv  
(4dhuT  
p=pr.'*pt; TwVlg;  
cM$P`{QrM  
%green's function 8>WC5%f*  
dAkgR~  
G=1; @jsDq Ln  
cs=1; T`;%TO*Y  
for i=1:1:m 8(~K~q[Cr  
  for j=1:1:n zhpt%7So  
    u(j,i)=heaviside(cs*t(i)-r(j))/(pi*G*sqrt(t(i)^2-(r(j)/cs)^2)); `m!j$,c.  
  end k=4N.*#`y  
end ~1W x=  
}}>q2y  
%convolution and response of displacement 32/MkuY^u  
,z-}t& _t  
for i=1:1:m K%F,='P}  
  for j=1:1:n AigS!-   
    v(j,i)=0; S/ODqL|  
    for k=1:1:i nysUZB  
      for g=1:1:j w6{TE(]zp  
        v(j,i)=v(j,i)+p(g,k)*u(j-g+1,i-k+1)*dr*dt; P#XID 2;  
      end; O]1y0BOQ  
    end; e0T34x'
  
  end; 1>E<8&2[L  
end; ZRg;/sX]  
RkBb$q9F]  
%plot the response history 3P^sM1  
m6[0Kws&  
plot(t,v(n,); s1h/}  
xlabel('t/s'); [N#,K02mk  
ylabel('v/m'); D-4f >  
grid on; :]^P1sH[  
title('response of point A or C'); NT+?#0I  
QUQu^p  
5. 直接用MATLAB自带卷积分算的出解（如下图5）及相应代码： 7lBAxqr2  
7w/4QiI  
pnbIiyV  
79xx2  
clear; EodQ*{l  
clc; pXtXjb  
w &(|e <  
%loading history f=mZu1(FZ  
O^^C;U@U<1  
dt=0.01; hnc@  
ti=0.0008; 0^RXGN  
te=12; zBk'{[y9L  
t=ti:dt:te; }Hxd*S  
pt=(10*t./1.5).*(t>=0 & t<=1.5)+(20-10*t./1.5).*(t>1.5 & t<=3)+0.*(t>3); 4bn(zyP  
m=length(pt); HY%i`]4X  
$l|qk z  
%load distribution in space ,~OwLWi-|X  
kT'u1q$3Vo  
rp=0.1; 0QOBL'{7)  
ri=0; @)4]b+8Z  
rc=0.1;  s8rE$  
dr=0.001; S7L=#+Z  
r=ri:dr:rc; Ksy -e{n  
pr=1*(r<=rp)+0*(r>rp); ,Qnd3[2[  
n=length(pr);  oze&  
H s"HID  
%load function with respect to t and r )>`G  
kMt 8/E`  
p=pr.'*pt; < 
VSA  
jhg;%+KB  
%green function 6w(6}m.L^  
A?t%e  
G=1; x*nSHb  
cs=1; ,}))u0q+:  
for i=1:1:m 5yiK+-iTs  
  for j=1:1:n OSf}Q=BL  
    u(j,i)=heaviside(cs*t(i)-r(j))/(pi*G*sqrt(t(i)^2-(r(j)/cs)^2)); Y~(Md@!0S  
  end <c,u3cp  
end 0Pe>Es|^A#  
HAJ7m!P  
%convolution and response of displacement 8peDI7[|  
\DD0s8  
v=conv2(u,p); V` 1/SQX  
x"{'&J[hx  
%plot the response history 2h=!k|6  
\,7}mdQSv  
rr=r(n); pD01,5/  
[tpu,rpu]=meshgrid(ti:dt:2*te-dt,ri:dr:2*rc); _Gjk;|Sx<I  
[X,Y,Z]=meshgrid(linspace(min(tpu(),max(tpu()),linspace(min(rpu(),max(rpu()),linspace(min(v(),max(v())); 66I"=:  
V=Y; [
.S#rGYk  
h=contourslice(X,Y,Z,V,tpu,rpu,v,[0 0]+rr); S4h:|jLUF  
set(h,'edgecolor','k'); S&\L-@  
contourslice(X,Y,Z,V,tpu,rpu,v,[0 0]+rr); .b-f9
qc=  
xlabel('t/s'); 2m35R&  
ylabel('r/m'); tP2qK_\e=  
zlabel('v/m'); YA +E\  
axis([0 12 0 0.1 0 12e4]); s+EAB{w$  
view(0,0); E8n)}[k!0  
grid on; 9J>&29@us0  
title('displacement response history of point A or C');