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

用MATLAB自带卷积分函数计算弹性半空间在一三角脉冲荷载下自由表面的竖向位移，与文献正确结果相差10的5次方倍，但自己编写卷积分代码，反而得到文献结果。请大家帮我看下问题出在哪？谢谢！ hgmCRC  
@~e5<:|5#  
另外，为何MATLAB自带卷积分函数结果矩阵的维数会是被卷积两矩阵维数之和减1？ -=="<0c  
+vH4MwG$.&  
1. 问题描述，如下图1所示： J,hCvm  
\+etCo   
#WuBL_nZ~  
u, ff>/1  
2. 格林函数，如下图2所示： s7<AfaJPF  
0sqFF[i  
>z03{=sAN  
]]mJ']l  
3. 参考文献正确解，如下图3所示： sK{e*[I>W  
9x8fhAy}4  
5R-6ji  
sB</DS  
4. 本人直接用自编卷积分MATLAB代码算的出解（如下图4）及相应代码： XSDpRo  
'%qr.T %  
R
i{=]$  
r$1Qf}J3=  
clear; |>Vb9:q9Po  
clc; )4O
xY[2J  
{=WgzP  
%loading history yfSmDPh  
d$RIS+V  
dt=0.01; `A>@]d  
ti=0.0008; rsQtMtS2  
te=12; -"`=1l  
t=ti:dt:te; 3mgD(,(^  
pt=(10*t./1.5).*(t>=0 & t<=1.5)+(20-10*t./1.5).*(t>1.5 & t<=3)+0.*(t>3); -@s#uA h  
m=length(pt); 3<!7>]A  
n]9$:aLZ  
%load distribution in space Ey2^?  
)UR7i8]!0  
rp=0.1; VRMXtQ*1Dm  
ri=0; x4 yR8n(  
rc=0.1; pb}*\/s  
dr=0.001; $X6h|?3U,  
r=ri:dr:rc; |N2#I
tBbW  
pr=1*(r<=rp)+0*(r>rp); Za9qjBH   
n=length(pr); tYS06P^<  
vt8By@]:  
%load function with respect to t and r n[z+<VGwC  
Wc#24:OKe3  
p=pr.'*pt; +2{Lh7Ks  
6t$8M[0-U  
%green's function qna8|3
eP  
Nc`L;CP  
G=1; Y|n"dMrL  
cs=1; "[J^YKoF  
for i=1:1:m +rd+0 `}C  
  for j=1:1:n e= AKD#  
    u(j,i)=heaviside(cs*t(i)-r(j))/(pi*G*sqrt(t(i)^2-(r(j)/cs)^2)); yA
t^;  
  end oxs#866x  
end ? k/`  
@5FQX  
%convolution and response of displacement bw7@5=?;  
t# i#(H  
for i=1:1:m b;n[mk  
  for j=1:1:n J zl6eo[;  
    v(j,i)=0; T[gv0|+  
    for k=1:1:i ]DcFySyv  
      for g=1:1:j HtFDlvdy]  
        v(j,i)=v(j,i)+p(g,k)*u(j-g+1,i-k+1)*dr*dt; $Yq9P0Ya  
      end; aOp\91  
    end; wT@og|M  
  end; icgfB-1|i  
end; b9krOe*j  
_b 0&!l<  
%plot the response history 6Oq7#3]  
HfVZ~PP  
plot(t,v(n,); +%'(!A?*`  
xlabel('t/s'); Da|z"I x  
ylabel('v/m'); mt .sucT  
grid on; qm}@!z^  
title('response of point A or C'); d0D]Q  
^!d3=}:
0  
5. 直接用MATLAB自带卷积分算的出解（如下图5）及相应代码： vN:Ng  
>6T8^Nt  
s<Ziegmw|g  
+>,I1{u%&  
clear; m`XH
KRp  
clc; 7
dWS  
qPNR`%}Q  
%loading history R_C)  
_f83-':W6  
dt=0.01; ^('wy};  
ti=0.0008; %EH)&k  
te=12; XSR 4iu  
t=ti:dt:te; V0@=^Bls  
pt=(10*t./1.5).*(t>=0 & t<=1.5)+(20-10*t./1.5).*(t>1.5 & t<=3)+0.*(t>3); # d  
m=length(pt); Vr}'.\$  
l#o  ~W`  
%load distribution in space .A
|udZ,  
S[gx{Bxiw  
rp=0.1; 7#XzrT]  
ri=0; qGo.WZ$  
rc=0.1; IxU/?Zm  
dr=0.001; 1zv'.uu.,  
r=ri:dr:rc; :;}P*T*PU  
pr=1*(r<=rp)+0*(r>rp); ?}oFg#m-<L  
n=length(pr); `?]k{ l1R  
la!~\wpa  
%load function with respect to t and r dPlV>IM$z  
}vuO$j  
p=pr.'*pt; RZLq]8pM  
FrS]|=LJhX  
%green function vXs"Dst  
tmqOJ  
G=1; ?s01@f#  
cs=1; [,Gg^*umS  
for i=1:1:m #mdc[.  
  for j=1:1:n o!Zb0/AP)  
    u(j,i)=heaviside(cs*t(i)-r(j))/(pi*G*sqrt(t(i)^2-(r(j)/cs)^2)); K+eM   
  end js(pC@<q5  
end u`W2+S  
SUiOJ[5,  
%convolution and response of displacement ftb\0,-  
(`^1Y3&2  
v=conv2(u,p); 04ui`-c(  
Lbgi7|&  
%plot the response history Wr 4,YQM  
pK*TE5]  
rr=r(n); 1EK*g;H  
[tpu,rpu]=meshgrid(ti:dt:2*te-dt,ri:dr:2*rc); dO'(2J8  
[X,Y,Z]=meshgrid(linspace(min(tpu(),max(tpu()),linspace(min(rpu(),max(rpu()),linspace(min(v(),max(v())); {: /}NpA$  
V=Y; 5m@V#2^P  
h=contourslice(X,Y,Z,V,tpu,rpu,v,[0 0]+rr); ?<!|  
set(h,'edgecolor','k'); oH@78D0A  
contourslice(X,Y,Z,V,tpu,rpu,v,[0 0]+rr); !$JT e  
xlabel('t/s'); #a#F,ZT  
ylabel('r/m'); KlEpzJ98  
zlabel('v/m'); O
-wzz  
axis([0 12 0 0.1 0 12e4]); x2xRBkRg=  
view(0,0); V3Bz Mw\9r  
grid on; Gc?a+T  
title('displacement response history of point A or C');