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

用MATLAB自带卷积分函数计算弹性半空间在一三角脉冲荷载下自由表面的竖向位移，与文献正确结果相差10的5次方倍，但自己编写卷积分代码，反而得到文献结果。请大家帮我看下问题出在哪？谢谢！ F[`dX  
Bx#=$ka  
另外，为何MATLAB自带卷积分函数结果矩阵的维数会是被卷积两矩阵维数之和减1？ %\n|2*r  
"Aw)0a[j1  
1. 问题描述，如下图1所示： A^A)arJS  
n${k^e-=  
bovAFdHW  
7mM
MVz2  
2. 格林函数，如下图2所示： .>P:{''  
cDE5/!  
Ym!e}`A\F  
qMA-#  
3. 参考文献正确解，如下图3所示： zNdkwj p+  
cC+2%q B  
4v3gpLH  

Pd(_  
4. 本人直接用自编卷积分MATLAB代码算的出解（如下图4）及相应代码： F"O\uo:3  
i. (Af$  
Ki7t?4YE  
<c:H u{D  
clear; 8N?D1;F;  
clc; slUi)@b  
h:r?:C>n  
%loading history SgehOu  
:Jv5Flxl  
dt=0.01; k+w J
i  
ti=0.0008; W I MBwmg  
te=12; u*rP8GuS  
t=ti:dt:te;
rDa{Ve  
pt=(10*t./1.5).*(t>=0 & t<=1.5)+(20-10*t./1.5).*(t>1.5 & t<=3)+0.*(t>3); 1 <+aF,  
m=length(pt); &>E gKL  
hqmE]hwc  
%load distribution in space j%^4 1y  
L/`1K_\l  
rp=0.1; x&0kIF'lq  
ri=0; Hq 3V+$  
rc=0.1; ]
sk=V.GGQ  
dr=0.001; |5O>7~Tp  
r=ri:dr:rc; ?+.C@_QZQ  
pr=1*(r<=rp)+0*(r>rp); +F2OPIanT~  
n=length(pr); lw.[qP  
2A[hMbL  
%load function with respect to t and r q CYu@Ho  
) ba~7A  
p=pr.'*pt; ?+^p$'5
 
1gbFl/i6T  
%green's function zyUS$g]&  
L\E>5G;  
G=1; l:uQ#Z)  
cs=1; IDFzyg_
 
for i=1:1:m +>K&zS  
  for j=1:1:n c@3 5\!9  
    u(j,i)=heaviside(cs*t(i)-r(j))/(pi*G*sqrt(t(i)^2-(r(j)/cs)^2)); Qz#By V:  
  end 7bihP@I!  
end Ve&_NVPrd  
f:<BUqa  
%convolution and response of displacement vZ"gCf3#?3  
FiUwy/,ZV  
for i=1:1:m "QxULiw  
  for j=1:1:n j-W$)c3X  
    v(j,i)=0; /#H P;>!n  
    for k=1:1:i ^jwzCo-  
      for g=1:1:j #?jsC)  
        v(j,i)=v(j,i)+p(g,k)*u(j-g+1,i-k+1)*dr*dt; ipbhjK$  
      end; #~"IlBk\  
    end; Y%;X7VxU*  
  end; fx[&"$X  
end; :TG;W,`.V  
Ez5
t)
l-  
%plot the response history zIjfxK  
J&,hC%]  
plot(t,v(n,); ~uty<fP  
xlabel('t/s'); XGH:'^o_  
ylabel('v/m'); kwc Cf2  
grid on; h-?yed*?  
title('response of point A or C'); Ohp@ZJ!a?  
TnK<Wba  
5. 直接用MATLAB自带卷积分算的出解（如下图5）及相应代码： 5~@-LXqL  
bS r"k  
xd^Pkf  
Y$A2{RjRq  
clear; 5P"R'/[PA_  
clc; iC=>wrqY>  
KGg S"d  
%loading history aSX4~UYB=  
X~0-WBz  
dt=0.01; Vb\g49\o/  
ti=0.0008; A#T"4'#?<  
te=12; 4^l9d  
t=ti:dt:te; n+ebi>
}P  
pt=(10*t./1.5).*(t>=0 & t<=1.5)+(20-10*t./1.5).*(t>1.5 & t<=3)+0.*(t>3); Pd"c*n&9  
m=length(pt); C1=&Vm>g+  
~io.TS|r  
%load distribution in space 8X"4RyNSn  
9$;5J  
rp=0.1; ~1wt=Ln>  
ri=0; tJrGRlB>  
rc=0.1; 0P9\;!Y  
dr=0.001; n-cI~Ax+4  
r=ri:dr:rc; fI<LxU
_n:  
pr=1*(r<=rp)+0*(r>rp); xw 43P.  
n=length(pr); az0=jou<Zl  
d\]KG(T  
%load function with respect to t and r v#%rjml[  
QhJN/v  
p=pr.'*pt; A+* lV*@0  
0lg'QG>  
%green function vjx'yh|  
+u0of^}=  
G=1; zdrP56rzZ  
cs=1; o?>0WSLlm  
for i=1:1:m 8:V,>PH  
  for j=1:1:n @tm2Y%Y!  
    u(j,i)=heaviside(cs*t(i)-r(j))/(pi*G*sqrt(t(i)^2-(r(j)/cs)^2)); z}u`45W+  
  end s)r!3HS  
end !~~KM?g  
`D~oY=  
%convolution and response of displacement  &kmaKc  
%"A8Af**I  
v=conv2(u,p); /-[vC$B"  
y 2> 93m  
%plot the response history p7;K] AW  
-@"3`uv"  
rr=r(n); PKrG6% W+  
[tpu,rpu]=meshgrid(ti:dt:2*te-dt,ri:dr:2*rc); 9d#?,:JG  
[X,Y,Z]=meshgrid(linspace(min(tpu(),max(tpu()),linspace(min(rpu(),max(rpu()),linspace(min(v(),max(v())); *pa hZiO  
V=Y; '`k7l7I[@  
h=contourslice(X,Y,Z,V,tpu,rpu,v,[0 0]+rr); yV:8>9wE8  
set(h,'edgecolor','k'); B]G2P`sN  
contourslice(X,Y,Z,V,tpu,rpu,v,[0 0]+rr); ue@/o,C>  
xlabel('t/s'); HJ7A/XW  
ylabel('r/m'); zMbFh_dcq  
zlabel('v/m'); 
M1-tRF  
axis([0 12 0 0.1 0 12e4]); qm!oJL  
view(0,0); ="& GU%$  
grid on; ;7:} iKU  
title('displacement response history of point A or C');