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

用MATLAB自带卷积分函数计算弹性半空间在一三角脉冲荷载下自由表面的竖向位移，与文献正确结果相差10的5次方倍，但自己编写卷积分代码，反而得到文献结果。请大家帮我看下问题出在哪？谢谢！ 5jQP"^g  
Fdw[CYHz  
另外，为何MATLAB自带卷积分函数结果矩阵的维数会是被卷积两矩阵维数之和减1？ 55
DzBV
  
/R
C!Yi  
1. 问题描述，如下图1所示： $ddYH  
:U q]~e  
I3Lsj}69  
_e_%U<\4  
2. 格林函数，如下图2所示： h%s  
w'0M>2  
Bg
h$P  
Ltw7b  
3. 参考文献正确解，如下图3所示： 0q>lW &J  
<`3(i\-X  
l{U3;  
@qDrTH]5  
4. 本人直接用自编卷积分MATLAB代码算的出解（如下图4）及相应代码： ,m?D\Pru  
@,&m`qzd+  
b1u'ukDP\  
{?y7'  
clear; E?mp6R]}%  
clc; xW9 s[X  
Q7
5^7Ga_  
%loading history X
gKG\C=3  
weV#%6=5\  
dt=0.01; +I n"OR%  
ti=0.0008; ewG21 q$  
te=12; 55LF  
t=ti:dt:te; \Ji2uGT  
pt=(10*t./1.5).*(t>=0 & t<=1.5)+(20-10*t./1.5).*(t>1.5 & t<=3)+0.*(t>3); \,!q[nC  
m=length(pt); +-\9'Q  
fti|3c  
%load distribution in space P` F'Nf2U  
I0vnd7  
rp=0.1; KaE;4gwM  
ri=0; C<t>m_t9  
rc=0.1; bW^QH-t  
dr=0.001; m#$za7  
r=ri:dr:rc; )JQQ4D  
pr=1*(r<=rp)+0*(r>rp); sri#L+I  
n=length(pr); F\R}no5C  
#6jwCEo=V  
%load function with respect to t and r cOZ^huK  
ebe@.ZVSi  
p=pr.'*pt; 1$VI\}  
_ICDtG^  
%green's function uW~,H}E  
j~H`*R=ld#  
G=1; x2sOEkcQ  
cs=1; E`n`#=xKR  
for i=1:1:m UMwMXmZNJ  
  for j=1:1:n ;cn.s,  
    u(j,i)=heaviside(cs*t(i)-r(j))/(pi*G*sqrt(t(i)^2-(r(j)/cs)^2)); ~ p.W*skD  
  end $jm<' 4  
end "3|"rc&F#  
Y^Q|l%Qrb  
%convolution and response of displacement 
bMZn7c  
X9A[  
for i=1:1:m g<4M!gi  
  for j=1:1:n |a$w;s>\  
    v(j,i)=0; 9sjW  
    for k=1:1:i <57l|}8  
      for g=1:1:j .GN$H>
')  
        v(j,i)=v(j,i)+p(g,k)*u(j-g+1,i-k+1)*dr*dt; 'F?Znd2L  
      end; "EYjY->  
    end;
!s*''v*  
  end; %`'z^W  
end; 8ysK VF  
)xx/di  
%plot the response history eJGos!>*  
K=?F3tX^  
plot(t,v(n,); &jZ|@K?  
xlabel('t/s'); W+ '}O<  
ylabel('v/m'); $cK B+}  
grid on; >Mz|e(6  
title('response of point A or C'); }!<cph  
w00\1'-Kz  
5. 直接用MATLAB自带卷积分算的出解（如下图5）及相应代码： _`{{39 F  
F` 5/9?;|  
5b`xN!c  
llfiNEK5;  
clear; 25c!-.5D  
clc; Z_ gVYa  
-xu.=n@,  
%loading history xO-U]%oq  
R(83E B~_  
dt=0.01; +7<>x-+  
ti=0.0008; ;T
{/;  
te=12; <lmJa#  
t=ti:dt:te; /)?P>!#;\  
pt=(10*t./1.5).*(t>=0 & t<=1.5)+(20-10*t./1.5).*(t>1.5 & t<=3)+0.*(t>3); So*Wk "  
m=length(pt); ynbpewaa  
,(27p6!  
%load distribution in space P&3/nL$9N  
~!-8l&C  
rp=0.1; N8YBu/  
ri=0; >DUE8hp;<  
rc=0.1; j~S!!Z]  
dr=0.001; K}<!{/fi)  
r=ri:dr:rc; fEG3b#t N  
pr=1*(r<=rp)+0*(r>rp); ,."(Gp  
n=length(pr); Gi2ad+QH-  
nl9Cdi]o  
%load function with respect to t and r Y0yO`W4  
:KP'xf.  
p=pr.'*pt; \seG2vw$  
AJ
` v  
%green function Rfc&OV  
AV 5\W}  
G=1; ]|t.wr3AU  
cs=1; )R)$T'  
for i=1:1:m E:4P1,%01+  
  for j=1:1:n 1R%`i'$/  
    u(j,i)=heaviside(cs*t(i)-r(j))/(pi*G*sqrt(t(i)^2-(r(j)/cs)^2)); u%
1k  
  end zO5u{  
end 8C,utjy  
$
%%>n^??  
%convolution and response of displacement BH0#Q5  

hAr[atu87  
v=conv2(u,p); LL[#b2CKa  
!8@rK$DB  
%plot the response history iynS4
]`U  
E}' d,v#Z{  
rr=r(n); EKd3$(^  
[tpu,rpu]=meshgrid(ti:dt:2*te-dt,ri:dr:2*rc); RvSq KW8  
[X,Y,Z]=meshgrid(linspace(min(tpu(),max(tpu()),linspace(min(rpu(),max(rpu()),linspace(min(v(),max(v())); rKQASRF5*  
V=Y; sMS9!{A  
h=contourslice(X,Y,Z,V,tpu,rpu,v,[0 0]+rr); px}7If  
set(h,'edgecolor','k'); V"by9p|V`  
contourslice(X,Y,Z,V,tpu,rpu,v,[0 0]+rr); $jed{N7Y  
xlabel('t/s'); hRa(<ZK  
ylabel('r/m'); 3).o"AN  
zlabel('v/m'); #f3;}1(  
axis([0 12 0 0.1 0 12e4]); 9X$#x90  
view(0,0); + lB+|yJ+  
grid on; bjPbl2K  
title('displacement response history of point A or C');