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

用MATLAB自带卷积分函数计算弹性半空间在一三角脉冲荷载下自由表面的竖向位移，与文献正确结果相差10的5次方倍，但自己编写卷积分代码，反而得到文献结果。请大家帮我看下问题出在哪？谢谢！ 75pn1*"gQ  
AJ>$`=  
另外，为何MATLAB自带卷积分函数结果矩阵的维数会是被卷积两矩阵维数之和减1？ 'l<Oj&E  
:-_"[:t 5Z  
1. 问题描述，如下图1所示： O7%8FY  
[!C!R$AMa  
|No9eZ8>.  
4p7j"d5  
2. 格林函数，如下图2所示： :IX,mDO  
DUSQh+C  
? o&goiM  
&_q8F,I \<  
3. 参考文献正确解，如下图3所示： (}5};v  
HVdB*QEH  
.HN4xL  
6i;q=N$'  
4. 本人直接用自编卷积分MATLAB代码算的出解（如下图4）及相应代码： t7yvd7  
LSR0yCU  
*~~J1.ja>  
iGSF5S  
clear; %ikPz~(  
clc; ]Exbuc  
KjMwrMgC  
%loading history U'^ G-@  
]Q
,&7D Ah  
dt=0.01; .}GOHW)}  
ti=0.0008; Pv=]7>e  
te=12; f9OY>|a9  
t=ti:dt:te; g[' 7$  
pt=(10*t./1.5).*(t>=0 & t<=1.5)+(20-10*t./1.5).*(t>1.5 & t<=3)+0.*(t>3); '0^lMQMg  
m=length(pt); Z`f?7/"B  
1g,Of
r  
%load distribution in space B}P!WRNmln  
fRxn,HyV  
rp=0.1; 7|"l/s9,  
ri=0; Y3#8]Z_"}O  
rc=0.1; W9{i~.zo  
dr=0.001; :]4s;q:m  
r=ri:dr:rc; IAWs}xIly  
pr=1*(r<=rp)+0*(r>rp); k&M~yb  
n=length(pr); \PD%=
~  
?VCp_Ji  
%load function with respect to t and r $> ;|  
/eT9W[a  
p=pr.'*pt; ]heVR&bQ  
vTo+jQs^  
%green's function bxPJ5oT  
OLWn0  
G=1; S(Z\h_m(  
cs=1; WL|71?@C  
for i=1:1:m q6hH]Q>w*  
  for j=1:1:n U# IPYyV  
    u(j,i)=heaviside(cs*t(i)-r(j))/(pi*G*sqrt(t(i)^2-(r(j)/cs)^2)); v-8{mK`9\  
  end ([|^3tM  
end LN)yQ-  
~c55LlO>  
%convolution and response of displacement ~Y{]yBGoF  
x[fp7*TiG  
for i=1:1:m 7L!}
F;yT  
  for j=1:1:n $WE=u9m  
    v(j,i)=0; roPC ^Q  
    for k=1:1:i PT~F^8,)  
      for g=1:1:j >Hmho'  
        v(j,i)=v(j,i)+p(g,k)*u(j-g+1,i-k+1)*dr*dt; me F.  
      end; y<~(}xsHh  
    end; X40
JCQx{+  
  end; 1;?w#/&t  
end; oGyo
U#z#  
}8ESp3~e_  
%plot the response history _+)n}Se  
4uH} SG[  
plot(t,v(n,); RameaFX8  
xlabel('t/s'); eoFG$X/PO  
ylabel('v/m'); dNCd-ep  
grid on; V6+:g=@U-l  
title('response of point A or C'); 4jlwu0L+  
BpGyjoJ2  
5. 直接用MATLAB自带卷积分算的出解（如下图5）及相应代码： p.<d+S<  
:?}>Q  
`9k\~D=D~  
,$aqF<+;  
clear; T24$lhM  
clc; 1
NG[   
F&#I[]#  
%loading history eL'fJcjw<  
Dw 5Ze  
dt=0.01; fB&i{_J  
ti=0.0008; zsj]WP6j  
te=12; z =\ENG|x#  
t=ti:dt:te; VZe'6?#  
pt=(10*t./1.5).*(t>=0 & t<=1.5)+(20-10*t./1.5).*(t>1.5 & t<=3)+0.*(t>3); DZ $O%  
m=length(pt); i+Mg[x$.  
|MN2v[y  
%load distribution in space qG2P?DR  
e|>@ >F]K  
rp=0.1; 9.,IqnP  
ri=0; 3g56[;Up?  
rc=0.1; RH$l?j6  
dr=0.001; *v: .]_;  
r=ri:dr:rc; 6ZwQ/~7H  
pr=1*(r<=rp)+0*(r>rp); nEP3B'+  
n=length(pr); _mQj
=  
DjiI*HLNR  
%load function with respect to t and r il"pKQF  
Z^
Wv(:Nr  
p=pr.'*pt; J9f]=1`  
. Y$xNLoP[  
%green function =EH/~NGk  
:T>OJ"p  
G=1; iA`.y9'2  
cs=1; arS@l<79  
for i=1:1:m 5E 9
R+N  
  for j=1:1:n X)=m4\R  
    u(j,i)=heaviside(cs*t(i)-r(j))/(pi*G*sqrt(t(i)^2-(r(j)/cs)^2)); c
/=\YeR  
  end n 4cos  
end HX,i{aWWy  
,1.Td=lY$  
%convolution and response of displacement Q'hs,t1<  
ooCfr?E  
v=conv2(u,p); +VJyGbOcC  
|9BX  ~`{  
%plot the response history ]GsI|se  
G)f!AuN=  
rr=r(n); Au{b1n  
[tpu,rpu]=meshgrid(ti:dt:2*te-dt,ri:dr:2*rc); F7/%,vf  
[X,Y,Z]=meshgrid(linspace(min(tpu(),max(tpu()),linspace(min(rpu(),max(rpu()),linspace(min(v(),max(v())); biForT_no  
V=Y; qP zxP @4  
h=contourslice(X,Y,Z,V,tpu,rpu,v,[0 0]+rr); z5D*UOy5M  
set(h,'edgecolor','k'); C[l5[DpH  
contourslice(X,Y,Z,V,tpu,rpu,v,[0 0]+rr); RE-y5.kE^  
xlabel('t/s'); sPl3JP&s  
ylabel('r/m'); gKmF#Z"\  
zlabel('v/m'); $Y\7E/T  
axis([0 12 0 0.1 0 12e4]); 3hEbM'L  
view(0,0); \/nSRAk  
grid on; ~]9EhC'l  
title('displacement response history of point A or C');