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

用MATLAB自带卷积分函数计算弹性半空间在一三角脉冲荷载下自由表面的竖向位移，与文献正确结果相差10的5次方倍，但自己编写卷积分代码，反而得到文献结果。请大家帮我看下问题出在哪？谢谢！ 3;/?q  
,+L KJl  
另外，为何MATLAB自带卷积分函数结果矩阵的维数会是被卷积两矩阵维数之和减1？ uDG+SdyN@  
>]$aoA#  
1. 问题描述，如下图1所示： SE`l(-tL  
(Pi-uL<[a  
(O5)wej   
YB!!/ SX4  
2. 格林函数，如下图2所示： >9(i)e  
(!zM\sF  
2_pz3<,\  
?b$3ob"  
3. 参考文献正确解，如下图3所示： :$H!@n*/R  
=Sxol>?t  
k$[{n'\@  
ZlR!s!vv  
4. 本人直接用自编卷积分MATLAB代码算的出解（如下图4）及相应代码： 'F_}xMU  
Aka^e\Y@6*  
cdp0!W4Gi  
1kFjas`g  
clear; D1"7s,Hmu  
clc; [8]m8=n  
RsSXhPk?  
%loading history gbGTG(:1S  
W"sr$K2m|  
dt=0.01; |O (G nsZ  
ti=0.0008; I-:`cON=G  
te=12; zXre~b03ZS  
t=ti:dt:te; it}-^3AM  
pt=(10*t./1.5).*(t>=0 & t<=1.5)+(20-10*t./1.5).*(t>1.5 & t<=3)+0.*(t>3); 1bRL"{m^)-  
m=length(pt); LpWI>sNv  
&4kM8Qh  
%load distribution in space j7/(sf  
#ooc)),  
rp=0.1; "bX4Q4Dq  
ri=0; f'{>AKi=C  
rc=0.1; !Yh}H<w0  
dr=0.001; kL7^$  
r=ri:dr:rc; kV)'a  
pr=1*(r<=rp)+0*(r>rp); HHS45kg[c  
n=length(pr); 'DAltr<  
*BOBH;s  
%load function with respect to t and r DX@}!6|T  
~mH+DV3  
p=pr.'*pt; FBYODw  
31XU7
A  
%green's function {+=i?  
olty4kGD$V  
G=1; `SOhG?Zo  
cs=1; S<oQ}+4[~  
for i=1:1:m {'~sS  
  for j=1:1:n iHz[Zw^.s  
    u(j,i)=heaviside(cs*t(i)-r(j))/(pi*G*sqrt(t(i)^2-(r(j)/cs)^2)); ,IjdO(?TC  
  end b=LF%P  
end \C/z%Hf7-  
M\UWWb&%\  
%convolution and response of displacement g_ M-F  
"{F;M{h$},  
for i=1:1:m ]h@{6N'oNS  
  for j=1:1:n njMLyT($  
    v(j,i)=0;  KOSyh<&  
    for k=1:1:i Q4%IxR?  
      for g=1:1:j p.Y$A if.  
        v(j,i)=v(j,i)+p(g,k)*u(j-g+1,i-k+1)*dr*dt; !<Z{@7oH  
      end; YvTA+yL  
    end; a$+#V=bA  
  end; 0j@IxEPs  
end; 8~5|KO >F  
9~Xg#{  
%plot the response history S}gD,7@  
;nk@XFJ  
plot(t,v(n,); 3?ba 1F0Nw  
xlabel('t/s'); |~NeB"l{  
ylabel('v/m'); .cR*P<3O
 
grid on; 2LhE]O(_"  
title('response of point A or C'); 79tJV  
QkX@QQT?  
5. 直接用MATLAB自带卷积分算的出解（如下图5）及相应代码： yiT{+;g^  
N$Hqa^!'T  
)BLmoJOf  
&&C~@WY,r  
clear; U42\.V0  
clc; eTZ`q_LfI1  
1g i}H)  
%loading history lIq~~cv)  
ay[+2"  
dt=0.01; O,9X8$5H-a  
ti=0.0008; +89o`u_l%  
te=12; |h,FUj<r  
t=ti:dt:te; 
N1? iiv  
pt=(10*t./1.5).*(t>=0 & t<=1.5)+(20-10*t./1.5).*(t>1.5 & t<=3)+0.*(t>3); oQvFrSz  
m=length(pt); AQ}l%  
v MWC(m  
%load distribution in space l<RfRqjw  
"k
>bUe|RG  
rp=0.1; \Da~p9T&  
ri=0; 6Bdyf(t  
rc=0.1; iEhDaC[e(b  
dr=0.001; b\L)m (  
r=ri:dr:rc; Yq;&F0paK  
pr=1*(r<=rp)+0*(r>rp); cEi<}9r  
n=length(pr); >B~?dTm  
a;p6?kv  
%load function with respect to t and r s1=u{ET  
dofR)"<p,^  
p=pr.'*pt; '3%*U*I  
Mf7E72{D  
%green function 7SHo%bA  
>sV Bj(f  
G=1; Gg+YfY_  
cs=1; Q-Y@)Mf~?0  
for i=1:1:m -A@U0=o  
  for j=1:1:n \UQ],+H  
    u(j,i)=heaviside(cs*t(i)-r(j))/(pi*G*sqrt(t(i)^2-(r(j)/cs)^2)); [+DNM 2A  
  end =g2\CIlVU6  
end ayH>XwY6  
)dg UmN
 
%convolution and response of displacement y''V"Be
  
\>[gl!B_Rr  
v=conv2(u,p); '%Dg{ zL  
M9g1d7%  
%plot the response history ZOHRUm  
@7|)RSBQz  
rr=r(n); yS"0/Rm}  
[tpu,rpu]=meshgrid(ti:dt:2*te-dt,ri:dr:2*rc); M,{<TpCx  
[X,Y,Z]=meshgrid(linspace(min(tpu(),max(tpu()),linspace(min(rpu(),max(rpu()),linspace(min(v(),max(v())); +~:0Dxv W  
V=Y; 6QptKXu7  
h=contourslice(X,Y,Z,V,tpu,rpu,v,[0 0]+rr); N7B}O*;  
set(h,'edgecolor','k'); s=jO;K$  
contourslice(X,Y,Z,V,tpu,rpu,v,[0 0]+rr); mt.,4  
xlabel('t/s'); uN&M\(  
ylabel('r/m'); WFdem/\kX  
zlabel('v/m'); v<fWc9
71  
axis([0 12 0 0.1 0 12e4]); 4H9xO[iM  
view(0,0); 2V<# Y  
grid on; Kz^hQd  
title('displacement response history of point A or C');