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

用MATLAB自带卷积分函数计算弹性半空间在一三角脉冲荷载下自由表面的竖向位移，与文献正确结果相差10的5次方倍，但自己编写卷积分代码，反而得到文献结果。请大家帮我看下问题出在哪？谢谢！ :S
ma`U&  
"]dI1 g_  
另外，为何MATLAB自带卷积分函数结果矩阵的维数会是被卷积两矩阵维数之和减1？ a =QCp4^  
z:;CX@)*  
1. 问题描述，如下图1所示： ,s(,
S  
HP=+<]?{G  
8_8l.!~  
=Uh$&m  
2. 格林函数，如下图2所示： xA
/D'  
RpF&\x>  
Ned
."e
  
KSvE~h[#+  
3. 参考文献正确解，如下图3所示： ys~x$  
7Wno':w8  
pUTr!fR  
OCUr{Nh  
4. 本人直接用自编卷积分MATLAB代码算的出解（如下图4）及相应代码： kl`W\tF  
HhpDR  
G?ZXWu.  
PdCEUh\>y  
clear; 6pzSp  
clc; (?c-iKGc  
/\Ef%@  
%loading history Fp:'MX  
q@qsp&0/  
dt=0.01; 99S^f:t  
ti=0.0008; eJSxn1GW  
te=12; g ?k=^C  
t=ti:dt:te; +H.`MZ=  
pt=(10*t./1.5).*(t>=0 & t<=1.5)+(20-10*t./1.5).*(t>1.5 & t<=3)+0.*(t>3); ]A"h&`Cvt  
m=length(pt); xmG<]WF>E  
T|p"0b A  
%load distribution in space .h[:xYm  
NgwbQ7)  
rp=0.1;
WM{=CD  
ri=0; R@0R`Zs  
rc=0.1; (=$x.1  
dr=0.001; g*Ph
v|kI  
r=ri:dr:rc; Q@niNDaW2  
pr=1*(r<=rp)+0*(r>rp); ^"g~-  
n=length(pr); B6"0OIDY"  
OPi

0~s  
%load function with respect to t and r /,dz@  
,>M[@4`,U  
p=pr.'*pt; 8QK&
_n*  
TJN4k@\$2  
%green's function P3%5?.S  
<[phnU^ 8  
G=1; ?CZdOl  
cs=1; sS Mh`4'  
for i=1:1:m H[gWGbPq7  
  for j=1:1:n GmG5[?)  
    u(j,i)=heaviside(cs*t(i)-r(j))/(pi*G*sqrt(t(i)^2-(r(j)/cs)^2)); [}:$yg  
  end U(Zq= M  
end Nl/dX-I  
9z0p5)]n>  
%convolution and response of displacement )7d&NE_  
phK/  
for i=1:1:m j [a(#V{  
  for j=1:1:n d1
*<Ll9K  
    v(j,i)=0; iwq!w6+  
    for k=1:1:i pJ=#zsE0  
      for g=1:1:j C}X\|J  
        v(j,i)=v(j,i)+p(g,k)*u(j-g+1,i-k+1)*dr*dt; }#J/fa9 !  
      end; #QPjkR|\  
    end; XuTD\g3)  
  end; qLCR] _*  
end; O8o3O 6[Y  
N;d] 14|  
%plot the response history
dI2 V>vk  
u y+pP!<  
plot(t,v(n,); y9;Yivr)  
xlabel('t/s'); /{[o~:'p  
ylabel('v/m'); ~dSr5LUD  
grid on; mR~&)QBP.  
title('response of point A or C'); ZG:{[sT  
: +u]S2u{  
5. 直接用MATLAB自带卷积分算的出解（如下图5）及相应代码： .6> w'F{>  
&L:!VL{I  
92c HwWZ!  
GVz6-T~\>  
clear; T+$[eWk"a  
clc; FlQGgVN  
~[ F`"  
%loading history @c#(.=  
)1z@
  
dt=0.01; >usL*b0%  
ti=0.0008; \GBuWY3B  
te=12; =v\.h=~~  
t=ti:dt:te; 43w}qY1  
pt=(10*t./1.5).*(t>=0 & t<=1.5)+(20-10*t./1.5).*(t>1.5 & t<=3)+0.*(t>3); LscGTs,  
m=length(pt); ,_P-$lB  
GB^Br6  
%load distribution in space b'y%n  
9$Y=orpWxr  
rp=0.1; lFkR=!?=  
ri=0; 9!GM{  
rc=0.1; 0%B/,/PxD  
dr=0.001; .d*8C,  
r=ri:dr:rc; CAlCDfKW}  
pr=1*(r<=rp)+0*(r>rp); HsWk*L `y  
n=length(pr); @d_M@\r=j  
QWU[@2@%r  
%load function with respect to t and r KXrjqqXs  
$:6!H:ty  
p=pr.'*pt; i@q&5;%%  
u?"Vm  
%green function )_:NLo:  
>ef6{UR
y<  
G=1; =%7-ZH9  
cs=1; *hrvYil2b  
for i=1:1:m Q
/?$x*\>  
  for j=1:1:n teP<!RKNb  
    u(j,i)=heaviside(cs*t(i)-r(j))/(pi*G*sqrt(t(i)^2-(r(j)/cs)^2)); [KQi.u  
  end ^pS~Z~[d/  
end Kq!3wb;  
3'Rx=G'  
%convolution and response of displacement }b}m3i1  
I'Hf{Erw  
v=conv2(u,p); jCY%|  
gr{ DWCK  
%plot the response history :]"V-1#}  
uNyVf7u  
rr=r(n); b$7 +;I;  
[tpu,rpu]=meshgrid(ti:dt:2*te-dt,ri:dr:2*rc); IgzQr>  
[X,Y,Z]=meshgrid(linspace(min(tpu(),max(tpu()),linspace(min(rpu(),max(rpu()),linspace(min(v(),max(v())); [WJ+h~~ o  
V=Y; FwK]$4*  
h=contourslice(X,Y,Z,V,tpu,rpu,v,[0 0]+rr); xLE)/}y_7H  
set(h,'edgecolor','k'); ,+VGSd  
contourslice(X,Y,Z,V,tpu,rpu,v,[0 0]+rr); 7^Uv7<pw  
xlabel('t/s'); h;Qk@F  
ylabel('r/m'); yu|>t4#GT  
zlabel('v/m'); d5b% W3  
axis([0 12 0 0.1 0 12e4]); dQvcXl]  
view(0,0); QPx^_jA  
grid on; t-AmX)$  
title('displacement response history of point A or C');