trl-008 发表于 2007-5-19 15:46

帮忙看看我二重积分的程序！！

function =jiexi(m,EI,K,C,E,I,Cs,P,r1,r2,v)
m=3.53*10^5;
EI=1.38*10^6;
K=1.7*10^8;
C=10^6;
E=1.5*10^9;
I=9.2*10^(-4);
Cs=6*10^9;
v=16.67;
P=4.9*10^7;
n=(Cs*I)/(2*m);
e=C/(2*m);
a=2*n*e-EI/m;
b=e^2-K/m;
r1=-a/(2*n^2)-1/n*sqrt((a/(2*n))^2-b);
r2=-a/(2*n^2)+1/n*sqrt((a/(2*n))^2-b);

if r1<0&r2>0
    syms w hx t;
    v1=-r2^(1/4);
    v2=r2^(1/4);
   
    for u=16.67:166.7:1667
   for v=1:10:100
      
    x=u
    t=v
    N=100;
    Q1=vectorize(char(exp(-(n.*w.^4+e).*(t-h)).*sin(sqrt(1./m.*(EI.*w.^4+K)-(n.*w.^4+e).^2).*(t-h)).*cos(w.*(x-v.*h))./sqrt(1./m.*(EI.*w.^4+K)-(n.*w.^4+e).^2)))
    Q2=vectorize(char(exp(-(n*w^4+e)*(t-h))*sin(sqrt((n*w^4+e)^2-1/m*(EI*w^4+K))*(t-h))*cos(w*(x-v*h))/sqrt((n*w^4+e)^2-1/m*(EI*w^4+K))))
    A1=dblquad(inline(Q1),0,v2,0,v)
    B1=dblquad(inline(Q2),v2,N,0,v)
    G1=A1+B1
    f=P/(pi*m)*G1
    end
f
end
f
谢谢啦，自己整了好久了，还是没弄出来，望高手们指点

trl-008 发表于 2007-5-19 16:14

回复 #1 trl-008 的帖子

高手们！！帮帮忙了
其中在第一次循环后是
A1 =

-2.2826e-004
第二次循环后
A1 =

      NaN +    NaNi
我单独拿出来编程积是具体的数，怎么回事？
而B1 =

      NaN +    NaNi

xjzuo 发表于 2007-5-19 16:51

请将问题用word贴一下．

trl-008 发表于 2007-5-19 17:06

谢谢高手！

xjzuo 发表于 2007-5-19 17:20

明显的一个问题：你的程序中显然应该引入i,j指标，以存储矩阵元．
代码也有些乱，我有时间再调试一下看看．

trl-008 发表于 2007-5-19 17:50

谢谢！！现在最大的问题积分中有问题！希望你帮我看看，小女子感激不尽！！:@P

trl-008 发表于 2007-5-20 09:03

回复 #5 xjzuo 的帖子

%移动载荷作用下粘弹性道路动力响应解析解
%给定mEI K C E I Cs=ET(T=0.05~5s) p的值
function =jiexi(m,EI,K,C,E,I,Cs,P,r1,r2,v)
m=3.53*10^5;
EI=1.38*10^6;
K=1.7*10^8;
C=10^6;
E=1.5*10^9;
I=9.2*10^(-4);
Cs=6*10^9;
v=16.67;
P=4.9*10^7;
n=(Cs*I)/(2*m);
e=C/(2*m);
a=2*n*e-EI/m;
b=e^2-K/m;
r1=-a/(2*n^2)-1/n*sqrt((a/(2*n))^2-b);
r2=-a/(2*n^2)+1/n*sqrt((a/(2*n))^2-b);
if r1<0&r2>0
    syms w hx t;
    v1=-r2^(1/4);
    v2=r2^(1/4);
    u=1:100:1001
v=1:10:101
for i=1:length(u)
      for j=1:length(v)
      
    x=u
    t=v
    N=100;
    Q1=vectorize(char(exp(-(n.*w.^4+e).*(t(j)-h)).*sin(sqrt(1./m.*(EI.*w.^4+K)-(n.*w.^4+e).^2).*(t(j)-h)).*cos(w.*(x(i)-v.*h))./sqrt(1./m.*(EI.*w.^4+K)-(n.*w.^4+e).^2)))
    Q2=vectorize(char(exp(-(n*w^4+e)*(t(j)-h))*sin(sqrt((n*w^4+e)^2-1/m*(EI*w^4+K))*(t(j)-h))*cos(w*(x(i)-v*h))/sqrt((n*w^4+e)^2-1/m*(EI*w^4+K))))
    A1(i,j)=dblquad(inline(Q1),0,v2,0,v)
    B1(i,j)=dblquad(inline(Q2),v2,N,0,v)
    G1(i,j)=A1(i,j)+B1(i,j)
    f(i,j)=P/(pi*m)*G1(i,j)
      end
end
surf(x,t,f)
hold on
grid on
xlabel('x')
ylabel('t')
zlabel('y')
title('函数图像')
end
高手，帮我看看啊，现在程序中加入了i,j指标，不过，直接第一个结果都出不来了:
:@L 急啊！

xjzuo 发表于 2007-5-20 09:30

那是因为你的写法明显有问题．
你的程序中有不少问题,我已经改正了.
问题是被积函数有奇异性,可能dblquad不适用，建议改为离散求和试试．
下面是我修改后正确的程序：
%%%======================================%%%
clear all
m=3.53*10^5;
EI=1.38*10^6;
K=1.7*10^8;
C=10^6;
E=1.5*10^9;
I=9.2*10^(-4);
Cs=6*10^9;
mu=16.67;　% 此处v字母被重复使用，已改正
P=4.9*10^7;
n=(Cs*I)/(2*m);
e=C/(2*m);
a=2*n*e-EI/m;
b=e^2-K/m;
r1=-a/(2*n^2)-1/n*sqrt((a/(2*n))^2-b);
r2=-a/(2*n^2)+1/n*sqrt((a/(2*n))^2-b);

if r1<0 & r2>0
    syms w h
    v1=-r2^(1/4);
    v2=r2^(1/4);
end　% 此处应该有一个 end

u=16.67:166.7:1667;
v=1:10:100;　% 此处v与上面的v字母重复

for i=1:length(u)
   for j=1:length(v)
    x=u(i);　
    t=v(j);　% 这部分修改后，注意下面的v也要改为t
    N=100;

Q1=exp(-(n.*w.^4+e).*(t-h)).*sin(sqrt(1./m.*(EI.*w.^4+K)-(n.*w.^4+e).^2).*(t-h)).*cos(w.*(x-mu.*h))./sqrt(1./m.*(EI.*w.^4+K)-(n.*w.^4+e).^2); 　% 什么一大堆的东西都可去掉

Q2=exp(-(n*w^4+e)*(t-h))*sin(sqrt((n*w^4+e)^2-1/m*(EI*w^4+K))*(t-h))*cos(w*(x-mu*h))/sqrt((n*w^4+e)^2-1/m*(EI*w^4+K));

    A1=dblquad(inline(Q1),0,v2,0,t);
    B1=dblquad(inline(Q2),v2,N,0,t);
    G1=A1+B1;
    f(i,j)=P/(pi*m)*G1;
end
end
surf(u,v,f)
%%%========================================%%%

[ 本帖最后由 xjzuo 于 2007-5-20 09:40 编辑 ]

trl-008 发表于 2007-5-20 14:57

回复 #8 xjzuo 的帖子

真的太厉害了，谢谢了

xiaoyuewei2001 发表于 2007-11-26 10:24

这个程序我在Matlab中怎么不能运行？？

??? Error using ==> inline.inline
Input must be a string.

Error in ==> jff at 37
    A1=dblquad(inline(Q1),0,v2,0,t);
页: [1]
查看完整版本: 帮忙看看我二重积分的程序！！