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

用MATLAB自带卷积分函数计算弹性半空间在一三角脉冲荷载下自由表面的竖向位移，与文献正确结果相差10的5次方倍，但自己编写卷积分代码，反而得到文献结果。请大家帮我看下问题出在哪？谢谢！ df\^uyD;  
i~Qnw-^B  
另外，为何MATLAB自带卷积分函数结果矩阵的维数会是被卷积两矩阵维数之和减1？ UHyGW$B  
qa-%j+  
1. 问题描述，如下图1所示： &t)$5\r  
jVlXB6[-  
&{4KymB:  
Q|KD$2rB  
2. 格林函数，如下图2所示： /]U),Lb
N  
{L'uuG\9U  
{
CH5`&  
%CoO-1@C  
3. 参考文献正确解，如下图3所示： )FQxVT,.  
z}BuR*WSY{  
F\
u]X  
M

r-l

 
4. 本人直接用自编卷积分MATLAB代码算的出解（如下图4）及相应代码： Vh?5  
GG
&J  
G*s5GG@Z.  
SI`ems{1>c  
clear; H0(.p'eN  
clc; E!A+J63zsw  
c1tM(]&  
%loading history >o:y.2yCe  
953G
mNZ7  
dt=0.01; vzX%x ul  
ti=0.0008; &s#OiF8  
te=12; |@W|
nbAfX  
t=ti:dt:te; J,G/L!Bp  
pt=(10*t./1.5).*(t>=0 & t<=1.5)+(20-10*t./1.5).*(t>1.5 & t<=3)+0.*(t>3); .R^R32ln  
m=length(pt); M{z&h>  
u+lNcyp"MW  
%load distribution in space /np05XhEa  
.(^%M 2:6  
rp=0.1; vRkVPkZ6|  
ri=0; ''^2
rF^  
rc=0.1; 73j\!x  
dr=0.001;
n +v(t  
r=ri:dr:rc; |zbM$37?k  
pr=1*(r<=rp)+0*(r>rp); a#D \8;  
n=length(pr);  sWyx_  
GvzaLEo  
%load function with respect to t and r 5Vc
~yMz  
$bU.6  
p=pr.'*pt; <=~*`eWV  
t/lQSUip  
%green's function l&|Tb8_'  
g es-nG-  
G=1; 8\F|{vt#  
cs=1; i);BTwW)#]  
for i=1:1:m `3eQ#,G!  
  for j=1:1:n #.<Dq8u  
    u(j,i)=heaviside(cs*t(i)-r(j))/(pi*G*sqrt(t(i)^2-(r(j)/cs)^2)); -G[TlH06  
  end \zh`z/=92  
end :]JMsa6  
Ts\PZQ!q  
%convolution and response of displacement vs^)=  
RD6>\9  
for i=1:1:m /H?) qk  
  for j=1:1:n 4`Cgz#v {  
    v(j,i)=0; zr ~4@JTS  
    for k=1:1:i '/s/o]'sUd  
      for g=1:1:j 5d;(D i5z  
        v(j,i)=v(j,i)+p(g,k)*u(j-g+1,i-k+1)*dr*dt; L)i6UAo  
      end; 9=J 3T66U  
    end; }#q0K  
  end; 8U
zF*gS  
end; Xz?7x0)Z  
+TW,!.NBG  
%plot the response history tUksIUYD\  
ba tXj]:  
plot(t,v(n,); >u\'k+=  
xlabel('t/s'); "r=p/"4D  
ylabel('v/m'); J8B0H1  
grid on; )j QrD`  
title('response of point A or C'); iu9+1+-  
,V9r2QY  
5. 直接用MATLAB自带卷积分算的出解（如下图5）及相应代码： HXN. ,[  
vA{DF{S4  
}tW1\@ =  
!gL1  
clear; G?^w <  
clc; z5_jx&^Z  
G%junS'zt  
%loading history usNq]  
TyvUdU  
dt=0.01; 2eRv{_  
ti=0.0008; 6>3zD)tG  
te=12; [x%8l,O #l  
t=ti:dt:te; eNK6=D|  
pt=(10*t./1.5).*(t>=0 & t<=1.5)+(20-10*t./1.5).*(t>1.5 & t<=3)+0.*(t>3); RA!8AS?  
m=length(pt); 4av  
)8taMC:H^  
%load distribution in space hltUf5m'b  
fo=@ X>S  
rp=0.1; pxI[/vS N  
ri=0; 6FX]b4  
rc=0.1; (tF/2cZk  
dr=0.001; |kF"p~s  
r=ri:dr:rc; T2A74>Nw  
pr=1*(r<=rp)+0*(r>rp); 8.&P4u i  
n=length(pr);
e< G[!m  
=eR#]d  
%load function with respect to t and r Ax 4R$P.]u  
~<}?pDA}~  
p=pr.'*pt; o{' JO3  
i)/#u+Y1P  
%green function (S?qxW?  
M<x><U#]A  
G=1; 5k c
?:U&  
cs=1; p m<K6I  
for i=1:1:m _ t.E_K  
  for j=1:1:n 4^*Z[6nt|  
    u(j,i)=heaviside(cs*t(i)-r(j))/(pi*G*sqrt(t(i)^2-(r(j)/cs)^2)); l$!Z};mw0E  
  end M=fhRCUB  
end ('`mPD,  
kaRjv   
%convolution and response of displacement l}FA&c"  
W6)XMl}n  
v=conv2(u,p); Bnz}:te}  
7H)tF
&  
%plot the response history ?IDkDv!na~  
x}f)P  
rr=r(n); _QtW)\)5\  
[tpu,rpu]=meshgrid(ti:dt:2*te-dt,ri:dr:2*rc); o9v.]tb  
[X,Y,Z]=meshgrid(linspace(min(tpu(),max(tpu()),linspace(min(rpu(),max(rpu()),linspace(min(v(),max(v())); !-7<x"avm  
V=Y; >J,IxRGi  
h=contourslice(X,Y,Z,V,tpu,rpu,v,[0 0]+rr); &m`@6\N(   
set(h,'edgecolor','k'); <899r \  
contourslice(X,Y,Z,V,tpu,rpu,v,[0 0]+rr); F)50 6
 
xlabel('t/s'); @sRb1+nn  
ylabel('r/m'); ?i\$U'2*z3  
zlabel('v/m'); t4p-pH'9b  
axis([0 12 0 0.1 0 12e4]); "/x/]Qx2  
view(0,0); rm<`H(cT  
grid on; Kww+lgzS  
title('displacement response history of point A or C');