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

用MATLAB自带卷积分函数计算弹性半空间在一三角脉冲荷载下自由表面的竖向位移，与文献正确结果相差10的5次方倍，但自己编写卷积分代码，反而得到文献结果。请大家帮我看下问题出在哪？谢谢！ c]SXcA;Pmv  
F7!g+LPc<  
另外，为何MATLAB自带卷积分函数结果矩阵的维数会是被卷积两矩阵维数之和减1？ {O ]^8#v^  
>/@wht4- j  
1. 问题描述，如下图1所示： $]v=2j  
V2As 5  
x3j)'`=15  
{t<E*5N]a  
2. 格林函数，如下图2所示： 3K2B7loD)~  
>yr:L{{D}G  
'MLp*3djF,  
l`S2bb6uMR  
3. 参考文献正确解，如下图3所示： F xek#  
@ev"{dY  

i2]7Bf)oV  
w;@NYMK)  
4. 本人直接用自编卷积分MATLAB代码算的出解（如下图4）及相应代码： h0m+u}oP_H  
#|<\q*<  
0_F6t-  
zl?Gd4  
clear; wKxw|Fpn  
clc; |\/`YRg>  
@d&(*9Y  
%loading history ]S@
zhQ  
[D
r'  
dt=0.01; <'n'>@  
ti=0.0008; BbU%p  
te=12; rC`pTN  
t=ti:dt:te; "DN0|%`M/  
pt=(10*t./1.5).*(t>=0 & t<=1.5)+(20-10*t./1.5).*(t>1.5 & t<=3)+0.*(t>3); }q<p;4<\F  
m=length(pt);
YfRjr  
BC:d@  
%load distribution in space &8p]yo2zO  
tB S+?N  
rp=0.1; '%Cc!63t*  
ri=0; /V-uo(n< .  
rc=0.1; <v&L90+s\;  
dr=0.001; H7{)"P]{f  
r=ri:dr:rc; O;zq(/,-l  
pr=1*(r<=rp)+0*(r>rp); 5`Q*  
n=length(pr); ,z4)A&F[c;  
lZ5LHUzP  
%load function with respect to t and r " pg5w  
'q3<R%^Q   
p=pr.'*pt; /e/%mo  
4eFqD;  
%green's function >A5*=@7bY?  
O52B  
G=1; JVfSmxy.  
cs=1; b|SDg%e  
for i=1:1:m srzlr-J  
  for j=1:1:n 'e&L53n  
    u(j,i)=heaviside(cs*t(i)-r(j))/(pi*G*sqrt(t(i)^2-(r(j)/cs)^2)); V{jQ=<)@e  
  end 0d>|2QV   
end ,7HlYPec  
.r ,wc*SF  
%convolution and response of displacement 7:[u.cd  
|7Dc7p"D  
for i=1:1:m voX4A pl  
  for j=1:1:n W&g@o@wa  
    v(j,i)=0; @:,B /B;  
    for k=1:1:i )RUx  
      for g=1:1:j #VM+.75o1  
        v(j,i)=v(j,i)+p(g,k)*u(j-g+1,i-k+1)*dr*dt; vIq>QXb;d  
      end; h^QLvOuR  
    end; A9 *P7  
  end; pc/x&VY%  
end; 8dPDs#Zl  
_HMQx_e0YM  
%plot the response history ?04jkq&  
xDqJsp=]-  
plot(t,v(n,); RSfB9)3D  
xlabel('t/s'); 91f{qq=#J{  
ylabel('v/m'); Nm{+!}cC  
grid on; a/wg%cWG_  
title('response of point A or C'); k|hy_? *  
V ,+&.A23  
5. 直接用MATLAB自带卷积分算的出解（如下图5）及相应代码： '6g;UOx^=  
=EJ8J;y_f  
53X i)  
h(~@ nd{  
clear; 1k)31GEQw  
clc; WswM5RN  
X%C`('"R  
%loading history ^X]rFY1  
*XUJv
&ZN  
dt=0.01; As{Q9o5j/  
ti=0.0008; uSsP'qd  
te=12; PF+F^;C  
t=ti:dt:te; ;*c8,I;  
pt=(10*t./1.5).*(t>=0 & t<=1.5)+(20-10*t./1.5).*(t>1.5 & t<=3)+0.*(t>3); 3VZ}5  
m=length(pt); %LM2CgH V  
h
5)4Z^n  
%load distribution in space b6%[?k  
$)YalZ  
rp=0.1; tDC?St1  
ri=0; r+Pfq[z&  
rc=0.1; =67ab_V  
dr=0.001; m|uVmg!*  
r=ri:dr:rc; FOyANN'  
pr=1*(r<=rp)+0*(r>rp); (Sv=R(_s  
n=length(pr); iv!;gMco  
p=XEMVqm  
%load function with respect to t and r Wq2Bo*[*  
"g
5<jp  
p=pr.'*pt; UuKW`(?^  

&gLXS1O  
%green function QbJE+m5  
7+NBcZuG9  
G=1; G1 K@Ir<  
cs=1; >b7Yk)[%  
for i=1:1:m c)j60y   
  for j=1:1:n uv|RpIve:  
    u(j,i)=heaviside(cs*t(i)-r(j))/(pi*G*sqrt(t(i)^2-(r(j)/cs)^2)); r@t \a+  
  end XQ'$J_hC  
end d53Eu`QW?  
/$\N_`bM  
%convolution and response of displacement Vg^yjP{sv  
-<6b[YA  
v=conv2(u,p); )H1\4LeP  
).32Im!;#R  
%plot the response history +xL' LCx  
>qy$W4  
rr=r(n);
Dz3~cuVb  
[tpu,rpu]=meshgrid(ti:dt:2*te-dt,ri:dr:2*rc); f8?hEa:js  
[X,Y,Z]=meshgrid(linspace(min(tpu(),max(tpu()),linspace(min(rpu(),max(rpu()),linspace(min(v(),max(v())); ]."t  
V=Y; SgWLs%B  
h=contourslice(X,Y,Z,V,tpu,rpu,v,[0 0]+rr); og?L 9  
set(h,'edgecolor','k'); 6ys|'<?  
contourslice(X,Y,Z,V,tpu,rpu,v,[0 0]+rr); <YM!K8hu$  
xlabel('t/s'); ]>/oo
=E  
ylabel('r/m'); ,.HS )<B  
zlabel('v/m'); l( WF  
axis([0 12 0 0.1 0 12e4]); J4]tT pu"K  
view(0,0); K_Pbzj4(P  
grid on; 5E#8F  
title('displacement response history of point A or C');