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

用MATLAB自带卷积分函数计算弹性半空间在一三角脉冲荷载下自由表面的竖向位移，与文献正确结果相差10的5次方倍，但自己编写卷积分代码，反而得到文献结果。请大家帮我看下问题出在哪？谢谢！ e 9p+  
X!'nfN  
另外，为何MATLAB自带卷积分函数结果矩阵的维数会是被卷积两矩阵维数之和减1？ Adyv>T9  
"~-Y'O  
1. 问题描述，如下图1所示： $d[ -feU  
e1d);m$  
qYi<GI*|@  
gr&Rkuyfv  
2. 格林函数，如下图2所示： <;T$?J9  
-( d,AX  
M?yWFqFt9m  
0SJ7Q
Ro|K  
3. 参考文献正确解，如下图3所示： 
CHZjK(a  
!"dn!X  
9[L@*7A`m  
Gb?O-z%8*  
4. 本人直接用自编卷积分MATLAB代码算的出解（如下图4）及相应代码： ww0m1FzX  
h3z=tu['  
1Vden.H*CI  
$ 2/T]  
clear; ,vN0Jpf}\8  
clc; i*q!
|^M  
Vv]81y15Q;  
%loading history 0lyCk}c  
HJV8P2f8`  
dt=0.01; qrq9NPf  
ti=0.0008; P2Or|_z  
te=12; ZJ|@^^GcL  
t=ti:dt:te; C/sDyv$  
pt=(10*t./1.5).*(t>=0 & t<=1.5)+(20-10*t./1.5).*(t>1.5 & t<=3)+0.*(t>3); ^KK9T5H  
m=length(pt); 8N58w)%7`  
HDTdOG)  
%load distribution in space m{ya%F  
Gkfc@[Z V  
rp=0.1; .W9/*cZV0  
ri=0; !edgziuO  
rc=0.1; Sn_zhQxG  
dr=0.001; tG{?  
r=ri:dr:rc; TLkJZ4}?Q  
pr=1*(r<=rp)+0*(r>rp); /
p&)bL  
n=length(pr); >Za66<:  
8G SO]R  
%load function with respect to t and r %5zztReI  
cv'Fc  
p=pr.'*pt; INHN=KY{  
0lvX,78G;  
%green's function VB?mr13}G  
H=7zd|W  
G=1; /_,} o7@t~  
cs=1; _z3Hl?qk=  
for i=1:1:m te+5@k#t  
  for j=1:1:n gUrb&#\X  
    u(j,i)=heaviside(cs*t(i)-r(j))/(pi*G*sqrt(t(i)^2-(r(j)/cs)^2)); a%wK[yVp  
  end #=MQE  
end h0N*hx   
d\cwUXf J  
%convolution and response of displacement K%p*:P  
Gn ]%'lrg'  
for i=1:1:m fGv`.T_d  
  for j=1:1:n ItoSORVV  
    v(j,i)=0; P'nby
F  
    for k=1:1:i 9t$%Tc#Z  
      for g=1:1:j GW(-'V/  
        v(j,i)=v(j,i)+p(g,k)*u(j-g+1,i-k+1)*dr*dt; Q)l]TgvSe  
      end; ^z[-pTY  
    end; (5"BKu1t  
  end; cZ" Ut  
end; $j~oB:3n7  
_n3Jf<Y  
%plot the response history AlQ!Q)y<@  
I:~L!%  
plot(t,v(n,); j=^b'dyL  
xlabel('t/s');
J6!t"eB+  
ylabel('v/m'); ;,z^!bD  
grid on; g>[|/z P  
title('response of point A or C'); +njE  
oadlyqlw#  
5. 直接用MATLAB自带卷积分算的出解（如下图5）及相应代码： n ^T_pqV?X  
TwZvz[u  
Yg;g!~
  
q5$z:'zE  
clear; %;.|?gR  
clc; %5_eos&<^)  
,u}n!quA  
%loading history EO|r  
zN\~v  
dt=0.01; NRS!Ox  
ti=0.0008; {C%/>e2-%  
te=12; _+,2b:D:  
t=ti:dt:te; `9QrkkG+  
pt=(10*t./1.5).*(t>=0 & t<=1.5)+(20-10*t./1.5).*(t>1.5 & t<=3)+0.*(t>3); dk
V%Pyj  
m=length(pt); !U"1ZsO)l  
(u]ajT  
%load distribution in space E(T6s^8  
TsPO+x$l  
rp=0.1; ;+-$=l3[a  
ri=0; ]|q\^k)JU  
rc=0.1; ,i2%FW  
dr=0.001; |Hbe]2"x>  
r=ri:dr:rc; ?l_>rSly5  
pr=1*(r<=rp)+0*(r>rp); mu1oD;lQ  
n=length(pr); pGi "*oZD  
;8~`fK  
%load function with respect to t and r XR^VRn6O  
vf@d(g  
p=pr.'*pt; 6e@ O88=  
AJrwl^lm  
%green function cU25]V^{\  
r\Wp\LfY&{  
G=1; j$*]'s&_hZ  
cs=1; XM/P2=;  
for i=1:1:m +a&-'`7g  
  for j=1:1:n ;G.m;5A
  
    u(j,i)=heaviside(cs*t(i)-r(j))/(pi*G*sqrt(t(i)^2-(r(j)/cs)^2)); g<s[6yA  
  end *@Z/L26s;=  
end ay2 m!s Q  
r'hr'wZ  
%convolution and response of displacement #R|M(Z">q  
`hM:U  
v=conv2(u,p); Ep}KIBBO  
O.=~/!(  
%plot the response history %E7+W{?*1  
:^SpKe(7  
rr=r(n); ->}K-n ),  
[tpu,rpu]=meshgrid(ti:dt:2*te-dt,ri:dr:2*rc); DYH-5yX7  
[X,Y,Z]=meshgrid(linspace(min(tpu(),max(tpu()),linspace(min(rpu(),max(rpu()),linspace(min(v(),max(v())); Z*kGWL  
V=Y; i:WHql"Kw_  
h=contourslice(X,Y,Z,V,tpu,rpu,v,[0 0]+rr); v@k62@;
 
set(h,'edgecolor','k'); ~?vm97l  
contourslice(X,Y,Z,V,tpu,rpu,v,[0 0]+rr); =JyYU*G4  
xlabel('t/s'); )2oWoZvi9  
ylabel('r/m'); 9`^VuC'  
zlabel('v/m'); Iz
2K  
axis([0 12 0 0.1 0 12e4]); 1!\!3xaV  
view(0,0); xIF z@9+k  
grid on; RlX;c!K  
title('displacement response history of point A or C');