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

用MATLAB自带卷积分函数计算弹性半空间在一三角脉冲荷载下自由表面的竖向位移，与文献正确结果相差10的5次方倍，但自己编写卷积分代码，反而得到文献结果。请大家帮我看下问题出在哪？谢谢！ h6(\ tRd!\  
5{Oq* |  
另外，为何MATLAB自带卷积分函数结果矩阵的维数会是被卷积两矩阵维数之和减1？ 7/969h^s  
N fBH  
1. 问题描述，如下图1所示： ?nCo?A  
u==`]\_@  
Pl\r|gS;  
q(9S4F   
2. 格林函数，如下图2所示： iU/v;T(  
iRIO~XVo  
2e<u/M21>  
6>Z)w}x^  
3. 参考文献正确解，如下图3所示： TCLXO0  
ecsQshR  
k E},>+W+  
Q
^{XM  
4. 本人直接用自编卷积分MATLAB代码算的出解（如下图4）及相应代码： >At* jg48  
'5r\o8RjN  
#7r13$>!  
8:sQB%BB  
clear; :?P>))vT%  
clc; C)?tf[!_6  
t}wwRWo2?f  
%loading history RA$%3L[A!  
) -^(Su(!  
dt=0.01; fWz=bJ"V  
ti=0.0008; 0Lx,qZ'  
te=12; +/n<]
?(T  
t=ti:dt:te; I R|[&}z  
pt=(10*t./1.5).*(t>=0 & t<=1.5)+(20-10*t./1.5).*(t>1.5 & t<=3)+0.*(t>3); #;])/8R%  
m=length(pt); xS+!/pBf"Y  
9q]n&5  
%load distribution in space rX}FhBl5  
{ u %xc"0y  
rp=0.1; :LxsiDrF[  
ri=0; )&s9QBo{b  
rc=0.1; %:!ILN  
dr=0.001;  &Sdf0"  
r=ri:dr:rc; 3%(,f,  
pr=1*(r<=rp)+0*(r>rp); ;-Ki`x.oJ  
n=length(pr); q
x1+'  
@:Emmzucv|  
%load function with respect to t and r ' +f(9/  
~$jRn(2  
p=pr.'*pt; r
cAPp  
Gq]/6igzX  
%green's function }Y!v"DO#Q*  
*_sSM+S  
G=1; 4Ifz-t/  
cs=1; Sfa;;7W@R  
for i=1:1:m \gFV6 H?`  
  for j=1:1:n nt_FqUJ  
    u(j,i)=heaviside(cs*t(i)-r(j))/(pi*G*sqrt(t(i)^2-(r(j)/cs)^2)); TZ#^AV=ae  
  end vyvb-oz;u  
end sH.,O9'r  
ix
_&os]L_  
%convolution and response of displacement Ke^9R-jP  
iyv
5\  
for i=1:1:m k![oJ.vHD  
  for j=1:1:n V<ii
 
    v(j,i)=0; hplxs#  
    for k=1:1:i O>eg_K,c  
      for g=1:1:j ,*.qa0E#W  
        v(j,i)=v(j,i)+p(g,k)*u(j-g+1,i-k+1)*dr*dt; t\WU}aKML  
      end; B8~bx%)3T  
    end; *M-'R*Np  
  end; 5c0$oyl)M  
end; wv  
(tCBbPW6T?  
%plot the response history >zfFvx_q  
&N*l?7(  
plot(t,v(n,); @:}la  
xlabel('t/s'); v,!`A!{D  
ylabel('v/m'); "0Z5cQjg  
grid on; &S39SV  
title('response of point A or C'); H6hhU'Kxf8  
#wZbG|%  
5. 直接用MATLAB自带卷积分算的出解（如下图5）及相应代码： E(_lm&,4+  
sD$K<nyz  
*z^Au7,&  
48_( 'z*>  
clear; hd W7Qck"  
clc; Dxe]LES\]  
b%].D(qBy  
%loading history u{cb[M  
bPIo9clq  
dt=0.01; K7t_Q8  
ti=0.0008; 0j2mTF(C  
te=12; j8+>E?nm  
t=ti:dt:te; 0jt@|
3  
pt=(10*t./1.5).*(t>=0 & t<=1.5)+(20-10*t./1.5).*(t>1.5 & t<=3)+0.*(t>3); hh[x(O)TC~  
m=length(pt); !p Q*m`Xo  
( 6ucA  
%load distribution in space LK<ZF=z]Z  
IEe;ygL#  
rp=0.1; OBf$Z"i  
ri=0; QT=i>X  
rc=0.1; 3G'cDemc  
dr=0.001; %I;uqf
 
r=ri:dr:rc; -EE}HUP)  
pr=1*(r<=rp)+0*(r>rp); >DAi
-`e  
n=length(pr); nG$+9}\UlP  
:1;"{=Yx}  
%load function with respect to t and r !AGoI7W}  
'
;m;K (g  
p=pr.'*pt; m95] z18T'  
L
b?0<  
%green function [OS&eK 8  
.hjN*4RY  
G=1; G[=;519  
cs=1; |d,bo/:  
for i=1:1:m .LGA0  
  for j=1:1:n %pLqX61t=  
    u(j,i)=heaviside(cs*t(i)-r(j))/(pi*G*sqrt(t(i)^2-(r(j)/cs)^2)); ,hZ?]P&  
  end ,BR W=  
end /4}y2JVv)  
So= BcX-  
%convolution and response of displacement 3_>=Cv}  
7d/I"?=|rA  
v=conv2(u,p); ANfy+@  
FH{p1_kZ=  
%plot the response history /]of@  
[
~kS)  
rr=r(n); ^O}J',Fm%f  
[tpu,rpu]=meshgrid(ti:dt:2*te-dt,ri:dr:2*rc); {"*_++|  
[X,Y,Z]=meshgrid(linspace(min(tpu(),max(tpu()),linspace(min(rpu(),max(rpu()),linspace(min(v(),max(v())); =}0$|@pl  
V=Y; Tfx-h)oP3  
h=contourslice(X,Y,Z,V,tpu,rpu,v,[0 0]+rr); 6n;? :./  
set(h,'edgecolor','k'); mC3:P5/c  
contourslice(X,Y,Z,V,tpu,rpu,v,[0 0]+rr); D~M*]&  
xlabel('t/s'); }*b\=AS=  
ylabel('r/m'); 8Un0<+b  
zlabel('v/m'); :(
4q\~  
axis([0 12 0 0.1 0 12e4]); )@<HG$#  
view(0,0); ~%h&ELSw  
grid on; 7X>*B~(R  
title('displacement response history of point A or C');