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

用MATLAB自带卷积分函数计算弹性半空间在一三角脉冲荷载下自由表面的竖向位移，与文献正确结果相差10的5次方倍，但自己编写卷积分代码，反而得到文献结果。请大家帮我看下问题出在哪？谢谢！ WnFG{S{s  
NIr@R7MKd  
另外，为何MATLAB自带卷积分函数结果矩阵的维数会是被卷积两矩阵维数之和减1？ 73A)lU.  
Bc-yxjsw  
1. 问题描述，如下图1所示： K[\'"HyQ,X  
UAF<m1  
}G46g#_6d>  
n@C~ev@%S  
2. 格林函数，如下图2所示： [36,eK  
W)j|rz.  
u]^N&2UW  
.Jb$l$5'w  
3. 参考文献正确解，如下图3所示： :yT-9Ze%q  
{)f~#37  
$5`!Z%>/  
ExSe=4q#  
4. 本人直接用自编卷积分MATLAB代码算的出解（如下图4）及相应代码： a\uie$"cr]  
5y2? f  
C8N{l:1f]  
iyZZ}
M  
clear; uNbH\qd=  
clc; ylf[/='0K  
C.:=lo B  
%loading history cR-~)UyrO  
U7mozHS,:9  
dt=0.01; ulHn#)  
ti=0.0008; _?7#MWe&  
te=12; ,''cNV  
t=ti:dt:te; C9n}6Er=,  
pt=(10*t./1.5).*(t>=0 & t<=1.5)+(20-10*t./1.5).*(t>1.5 & t<=3)+0.*(t>3); jg 2qGC  
m=length(pt); :A46~UA!$  
z!Q
DTIb  
%load distribution in space :^ i9]  
`+lHeLz':  
rp=0.1; XALI<ZY  
ri=0; p_*M:P1Ma4  
rc=0.1; *MNHT`Y^o  
dr=0.001; ~d{.ng 4K  
r=ri:dr:rc; =!Vf  
pr=1*(r<=rp)+0*(r>rp); M_0zC1  
n=length(pr); g o5]<4`r  
1xNVdI  
%load function with respect to t and r d&cU*  
:R6
bq!  
p=pr.'*pt; SQsSa1  
jcCoan  
%green's function QlFZO4 P3|  
\hO2p6  
G=1; +YOKA*  
cs=1; ?zJpD8e  
for i=1:1:m y<R=  
  for j=1:1:n rRES8/  
    u(j,i)=heaviside(cs*t(i)-r(j))/(pi*G*sqrt(t(i)^2-(r(j)/cs)^2)); PeX1wK%f  
  end 4W4kwU6D  
end &MR/6"/s  
@Fv=u  
%convolution and response of displacement z9 u$~  
){s*n=KIO  
for i=1:1:m /il@`w;G  
  for j=1:1:n
vqslirC  
    v(j,i)=0; #yseiVm;  
    for k=1:1:i !U_K&f  
      for g=1:1:j ;P&y,:<m:  
        v(j,i)=v(j,i)+p(g,k)*u(j-g+1,i-k+1)*dr*dt; sH,kW|D  
      end; 6TWWlU^e  
    end; 1 "TVRb  
  end; 5/[H+O1;  
end; =6FUNvP#8  
1PaUI#X"2F  
%plot the response history }y%`)lz~;  
A\rt6/  
plot(t,v(n,); :H6FPV78  
xlabel('t/s'); ,7Y-k'7Kop  
ylabel('v/m'); &WXY'A=  
grid on; a~h:qpgc  
title('response of point A or C'); E9j+o y  
z@s5m}  
5. 直接用MATLAB自带卷积分算的出解（如下图5）及相应代码： IJOvnZ("A  
O40+M)e]  
rn@`yTw^  
`"yxdlXA  
clear; JN/UUfj  
clc; y #f QPR  
?q`0ZuAg\<  
%loading history wo2@hav  
c;f!!3&  
dt=0.01; z_;3H,z`  
ti=0.0008; ymY1o$qWB}  
te=12; _eSdnHWx  
t=ti:dt:te; 5OIc(YhYf  
pt=(10*t./1.5).*(t>=0 & t<=1.5)+(20-10*t./1.5).*(t>1.5 & t<=3)+0.*(t>3); K)7zKEp`cj  
m=length(pt); q:>^ "P{  
n>,L=wV  
%load distribution in space o 6{\Zzp  
'PZ|:9FX!  
rp=0.1; Bsf7mcXz7z  
ri=0; 9DQ)cy  
rc=0.1; pN6%&@) =  
dr=0.001; >$677  
r=ri:dr:rc; x"kjs.d7[<  
pr=1*(r<=rp)+0*(r>rp); >t,M  
n=length(pr); D\~zS`}  
Gz I~TWc+G  
%load function with respect to t and r 14eW4~Mr  
+j+5ud`  
p=pr.'*pt; os3 8u!3-  
uxn)R#?  
%green function ?[TfpAtQ`  
d|9b~_::V  
G=1; dCYCHHHF  
cs=1; PW(\4Q\  
for i=1:1:m %OR|^M  
  for j=1:1:n 09KcKhFB  
    u(j,i)=heaviside(cs*t(i)-r(j))/(pi*G*sqrt(t(i)^2-(r(j)/cs)^2)); )CPM7>  
  end RhI;;
Y#@  
end CF!Sa6  
psh^MX)Q  
%convolution and response of displacement MmPU7Nl%X  
cxeghy:;U  
v=conv2(u,p); *F^wtH`  
3:/'t{ ^B  
%plot the response history 9L0GLmLk1u  
:6J +%(f  
rr=r(n); t2
2;87&|  
[tpu,rpu]=meshgrid(ti:dt:2*te-dt,ri:dr:2*rc); !9*c8bL D  
[X,Y,Z]=meshgrid(linspace(min(tpu(),max(tpu()),linspace(min(rpu(),max(rpu()),linspace(min(v(),max(v())); I:&/`K4,x,  
V=Y; A*h{Lsx;  
h=contourslice(X,Y,Z,V,tpu,rpu,v,[0 0]+rr); 3H\w2V  
set(h,'edgecolor','k'); *YTo{~  
contourslice(X,Y,Z,V,tpu,rpu,v,[0 0]+rr); h<<>
3A  
xlabel('t/s'); p1pQU={<  
ylabel('r/m'); kB:Uu}(=N  
zlabel('v/m'); @K223?c8l  
axis([0 12 0 0.1 0 12e4]); 1[F3 Z  
view(0,0); Y&H}xn  
grid on; _i_Q?w`  
title('displacement response history of point A or C');