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

用MATLAB自带卷积分函数计算弹性半空间在一三角脉冲荷载下自由表面的竖向位移，与文献正确结果相差10的5次方倍，但自己编写卷积分代码，反而得到文献结果。请大家帮我看下问题出在哪？谢谢！ yVp,)T9  
yM`u]p1  
另外，为何MATLAB自带卷积分函数结果矩阵的维数会是被卷积两矩阵维数之和减1？ Ta[\BWR2  

\_  
1. 问题描述，如下图1所示： )3)7zulnXH  
p="0Y<2l  
C5~n^I|  
X
o:Mar  
2. 格林函数，如下图2所示： `yP`
5a/  
xKl1DIN[  
Vef!5]t5  
M!kSt1  
3. 参考文献正确解，如下图3所示： )o<^6Ic%7  
@H<*|3J  
KIcIYCBz  
OLqV#i[K#9  
4. 本人直接用自编卷积分MATLAB代码算的出解（如下图4）及相应代码： tXqX[Td`0g  
b8_F2  
r ^\(M {  
|j-ng;  
clear; "X^<g{]  
clc; L
P=y$B  
T-#4hY`  
%loading history "g)V&Lx#X  
="5D}%  
dt=0.01; eXMIRus(  
ti=0.0008; c6
lCF &  
te=12; ZSs@9ej  
t=ti:dt:te; @({=~ W^  
pt=(10*t./1.5).*(t>=0 & t<=1.5)+(20-10*t./1.5).*(t>1.5 & t<=3)+0.*(t>3); $C sE[+k1  
m=length(pt); op\$(7<d-  
Rpv[rvK'  
%load distribution in space qY$ [2]  
0-[naGz  
rp=0.1; NYr)=&)Ke.  
ri=0; 4=ovm[  
rc=0.1; cS'{h  
dr=0.001; l

^d'8n
 
r=ri:dr:rc; qDG2rFu&[  
pr=1*(r<=rp)+0*(r>rp); >[Wjzg  
n=length(pr); B;?)X&n|X  
=@0J:"c  
%load function with respect to t and r + c+i u6+"  
Aj|->Y  
p=pr.'*pt; 2lO(f+  
)|vy}Jf7  
%green's function ^86M94k  
s[sv4hq  
G=1; f@$W5*j  
cs=1; 8]Tv1Wc  
for i=1:1:m l -xc*lC  
  for j=1:1:n i%JJ+9N  
    u(j,i)=heaviside(cs*t(i)-r(j))/(pi*G*sqrt(t(i)^2-(r(j)/cs)^2)); x1?mE)
n]  
  end Ix6\5}.c9  
end K"x_=^,Yu*  

pr,,E[  
%convolution and response of displacement <gFa@at  
"p*'H
Q  
for i=1:1:m lcm3wJ'w  
  for j=1:1:n tfN[-3)Z  
    v(j,i)=0; E*u*LMm  
    for k=1:1:i J8!2Tt  
      for g=1:1:j xnt)1Q  
        v(j,i)=v(j,i)+p(g,k)*u(j-g+1,i-k+1)*dr*dt; ?
].MnwYo  
      end; 9(J,&)J  
    end; 8}m J)9<7  
  end; e6>G8d  
end; &92/qRh7  
e`S\-t?Z  
%plot the response history +]nIr'V  
I8YUq  
plot(t,v(n,); DnFzCJ  
xlabel('t/s'); D;NL*4zt  
ylabel('v/m'); TIxOMYy  
grid on; ^?Mp(o  
title('response of point A or C'); ROmmak(y8  
*!ng)3#  
5. 直接用MATLAB自带卷积分算的出解（如下图5）及相应代码： 1C8xJ6F  
Ps>:|j+  
n."n?C'{  
l'/`2Y1  
clear; 0F@~[W|2  
clc; W2#
<]]-  
_ ,s^  
%loading history [#C6K '  
FGx)?  
dt=0.01; x};~8lGT>t  
ti=0.0008; 0U@#&pUc  
te=12; QM#Vl19>j(  
t=ti:dt:te; }L)[>  
pt=(10*t./1.5).*(t>=0 & t<=1.5)+(20-10*t./1.5).*(t>1.5 & t<=3)+0.*(t>3); 6e rY
jq  
m=length(pt); $3PDe  
]c~yMA+]FZ  
%load distribution in space `hZh}K^  
>7PQOQMW'
 
rp=0.1; 9xO@_pkX  
ri=0; MzX&|wimb  
rc=0.1; v@soS1V!  
dr=0.001; aZ4?!JW.  
r=ri:dr:rc; CNefk$/cR  
pr=1*(r<=rp)+0*(r>rp); 8W{M}>;[9  
n=length(pr); ^hiY6N &  
[B#XA}w  
%load function with respect to t and r Jz~:  
9zb1t1[W  
p=pr.'*pt; !9WGZfK+0Y  
$Yt|XT+!&  
%green function k9f|R*LM  
0M"n  
G=1; (0H=f6N  
cs=1; pX{wEc6}  
for i=1:1:m |y[I!JdR  
  for j=1:1:n *qm|A{FQR  
    u(j,i)=heaviside(cs*t(i)-r(j))/(pi*G*sqrt(t(i)^2-(r(j)/cs)^2)); NX,m6u  
  end CYLab5A  
end v>#Njgo  
P*A+k"DU1  
%convolution and response of displacement jkx>o?s)z  
Yu\$Y0 {]  
v=conv2(u,p); jel:oy|_  
Zs />_w}  
%plot the response history ?UAuUFueA  
YD'gyP4  
rr=r(n); dI ,A;.  
[tpu,rpu]=meshgrid(ti:dt:2*te-dt,ri:dr:2*rc); Ba@~:  
[X,Y,Z]=meshgrid(linspace(min(tpu(),max(tpu()),linspace(min(rpu(),max(rpu()),linspace(min(v(),max(v())); gns}%\,  
V=Y; aL4^ po  
h=contourslice(X,Y,Z,V,tpu,rpu,v,[0 0]+rr); (<r)xkn  
set(h,'edgecolor','k'); aWIkp5BFj  
contourslice(X,Y,Z,V,tpu,rpu,v,[0 0]+rr); &J
&'J~N  
xlabel('t/s'); T2MC`s|`  
ylabel('r/m'); o6px1C:  
zlabel('v/m'); I09 W=  
axis([0 12 0 0.1 0 12e4]); {#qUZ z-  
view(0,0); B`Z3e%g#  
grid on; zPa2fS8  
title('displacement response history of point A or C');