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

用MATLAB自带卷积分函数计算弹性半空间在一三角脉冲荷载下自由表面的竖向位移，与文献正确结果相差10的5次方倍，但自己编写卷积分代码，反而得到文献结果。请大家帮我看下问题出在哪？谢谢！ qG?Qc (  
-w}]fb2Q>  
另外，为何MATLAB自带卷积分函数结果矩阵的维数会是被卷积两矩阵维数之和减1？ G\Cp7:j}  
4C61GB?Vy  
1. 问题描述，如下图1所示： lhAX;s&9  
t(NI-UXBp  
t\~P:"  
g(qJN<RC/  
2. 格林函数，如下图2所示： Oj3.q#)`Z  
7vrl'^1  
{GK;63`1  
P3x= 8_#  
3. 参考文献正确解，如下图3所示： ff,pvk8N5  
75f"'nJ)  
"/3'XOK|  
diL+:H  
4. 本人直接用自编卷积分MATLAB代码算的出解（如下图4）及相应代码： @s ?  
[65`$x
-  
M/>7pZW  
-.u]Ge
My  
clear; t^R][Ay&  
clc; :t8b39
  
bnq;)>&  
%loading history (:TjoXXiY  
)NXmn95  
dt=0.01; F;4vPbH+  
ti=0.0008; tl,.fjZn  
te=12; .f%fHj  
t=ti:dt:te; k;AD`7(=  
pt=(10*t./1.5).*(t>=0 & t<=1.5)+(20-10*t./1.5).*(t>1.5 & t<=3)+0.*(t>3); Wz49
i9e+d  
m=length(pt); Sq/ qu-%X  
[q
)8N  
%load distribution in space &_dt>.  
-:Da&V  
rp=0.1; ; >hNt  
ri=0; 4:$4u@   
rc=0.1; (2J: #  
dr=0.001; QwJVS(Gs4  
r=ri:dr:rc; eg\v0Y!rI  
pr=1*(r<=rp)+0*(r>rp); Pq;U&,  
n=length(pr); 
aQ?/%\>  
la0BiLzb]  
%load function with respect to t and r "GMBjT8  
([T>.s  
p=pr.'*pt; P;=n9hgHI  
|:nOp(A\*  
%green's function O`x;,6Vr  
5cL83FQh  
G=1; q<[P6
}.  
cs=1; /YW>*?"N  
for i=1:1:m Wuc S:8#|  
  for j=1:1:n
8<S~Z:JK  
    u(j,i)=heaviside(cs*t(i)-r(j))/(pi*G*sqrt(t(i)^2-(r(j)/cs)^2)); ZM!CaR  
  end lYVz3p  
end oTU!R ,  
}Z@ovsG  
%convolution and response of displacement r3&G)g=u  
9ifDcYl  
for i=1:1:m |[<_GQl  
  for j=1:1:n *4Thd:7 `  
    v(j,i)=0; rb5~XnJk  
    for k=1:1:i =n5zM._S-  
      for g=1:1:j \o}xF@sM5  
        v(j,i)=v(j,i)+p(g,k)*u(j-g+1,i-k+1)*dr*dt; BP'36?=Zo  
      end; +%T\`6  
    end; vj'wm}/  
  end; qT{U(  
end; 3G,Oba[$<  
z&#SPH*  
%plot the response history [YF>:ydk  
8uc1iB  
plot(t,v(n,); `w#Oih!6A|  
xlabel('t/s'); R]c+?4J  
ylabel('v/m'); l&O
KBUG  
grid on; I5 o)_nc  
title('response of point A or C'); [842&5Pd?  
X$ 0?j1  
5. 直接用MATLAB自带卷积分算的出解（如下图5）及相应代码： VRWAm>u  
c}Ft^Il  
fi-WZ  
OE_XCZ!5P  
clear; a oD`=I*<  
clc; LSa,1{  
z1PBM
SG  
%loading history p4.wh|n  
]/[FR5>  
dt=0.01; jSh5!6O
  
ti=0.0008; m[?E  
te=12; ddJQC|xR}  
t=ti:dt:te; Vwg
|K|  
pt=(10*t./1.5).*(t>=0 & t<=1.5)+(20-10*t./1.5).*(t>1.5 & t<=3)+0.*(t>3); xu/cq9  
m=length(pt); h58`XH  
1an^1!  
%load distribution in space T! Y@`Ox  
,&]S(|2%>t  
rp=0.1; .zA^)qgL  
ri=0; H*RC@O_hv  
rc=0.1; twL3\ }N/B  
dr=0.001; 0%9 q8M;  
r=ri:dr:rc; BgurzS4-  
pr=1*(r<=rp)+0*(r>rp); >Wm`v.-  
n=length(pr); _E&A{HkJ  
q8X feoUV  
%load function with respect to t and r  8n#HFJ~  
Xb:;</  
p=pr.'*pt; :1cV;gJ  
c]x1HvPE  
%green function gn8R[5:!V  
<Swt);  
G=1; Uol|9F  
cs=1; Qi,j+xBp  
for i=1:1:m B:b5UD  
  for j=1:1:n Ygm`ZA y  
    u(j,i)=heaviside(cs*t(i)-r(j))/(pi*G*sqrt(t(i)^2-(r(j)/cs)^2)); \\r)Ue]  
  end eJF5n#  
end 2Nu=/tMN  
?i7%x,g(Z  
%convolution and response of displacement hm84Aq= f  
Y>|B;Kj0(  
v=conv2(u,p); t
X9{hC^  
|{BIHgMh  
%plot the response history *xx'@e|<;  
5gH1.7i b  
rr=r(n); X[*<NN  
[tpu,rpu]=meshgrid(ti:dt:2*te-dt,ri:dr:2*rc); 1tEgl\u\  
[X,Y,Z]=meshgrid(linspace(min(tpu(),max(tpu()),linspace(min(rpu(),max(rpu()),linspace(min(v(),max(v())); FOv=!'So  
V=Y; 8{wwd:6  
h=contourslice(X,Y,Z,V,tpu,rpu,v,[0 0]+rr); *W4m3Lq  
set(h,'edgecolor','k'); 9oRy)_5Z(=  
contourslice(X,Y,Z,V,tpu,rpu,v,[0 0]+rr); w k(VR  
xlabel('t/s'); lGV0*Cji  
ylabel('r/m'); _X^1IaL  
zlabel('v/m'); oX#Q<2z*  
axis([0 12 0 0.1 0 12e4]); ^=BTz9QM  
view(0,0); fM]+SMZy  
grid on; 63q^ $I  
title('displacement response history of point A or C');