应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 计算模拟 » 很简单的一个matlab程序运行时间过长，得不到结果
51/1返回列表
查看: 4935  |  回复: 4
只看楼主@他人 存档 新回复提醒 (忽略) 收藏    在APP中查看

hustboshao

铜虫 (小有名气)

[求助] 很简单的一个matlab程序运行时间过长，得不到结果已有1人参与

主程序：
clc
clear all
global kesai omega0 m rou D a1 a3 K C %%主系统基本参数
global epsilong_T s_33 C_pzt C_pzt1 K_pzt K_pzt1 k0 d_33 s_0 l_0 theta L R %压电材料参数
kesai=0.0013;m=0.44;omega0=62.83;rou=1.2;D=0.015;
a1=2.3;a3=-18;
K=1736.947916;
C=0.07187752;
s_0=0.005*0.005*pi; %cross section area of the PZT rod
l_0=0.01; %length of the PZT rod
epsilong_T=3400*8.854e-12; % permittivity under constant stress 常应力下的介电常数
s_33=20.7e-12; %compliance of the piezoceramic under constant electric 恒定电场下的弹性柔顺系数
d_33=593e-12; %压电常数
C_pzt=epsilong_T*s_0/l_0; %capacitance between the electrodes of the PZT rod with no external force无外力下的电容
K_pzt=(1/s_33)*s_0/l_0; %stiffness of the short-circuited PZT rod 短路刚度
k0=d_33*(1/(sqrt(epsilong_T*s_33))); % electromechanical coupling cofficient 机电耦合系数
%C_pzt1=C_pzt*(1-k0^2);%capacitance of the PZT rod under constant strain 常应变下的电容
C_pzt1=1.2e-7;
%K_pzt1=K_pzt*(1/(1-k0^2)); % stiffness of the PZT rod with open electrodes 回路刚度
%K_pzt1=C_pzt1/(epsilong_T*s_33*(1-k0^2)^2);
K_pzt1=K*0.25;
%theta=k0/(1-k0^2)*(1/(sqrt(epsilong_T*s_33)));% electromechanical coupling factor
theta=1.55e-3;
L=0.000001; % 电感
R=10; % 电阻
step_size=0.01;  %时间步长
precision=1e-9;   %精度量级1
options=odeset('reltol',precision);
N=20;
NN=N/0.01;
tic
U=5
%[T1,Y1]=ode45(@uncontrol,[0:step_size:N],[0.001 0],options,U);
[T2,Y2]=ode45(@control,[0:step_size:N],[0.001 0 0.001 0],options,U);
toc
%w_1=max(Y1(NN/2:end,1))/0.015
w_2=max(Y2(NN/2:end,1))
w_2_Q=max(Y2(NN/2:end,3))
%figure(1)
%plot(T1,Y1(:,1)/0.015,'r')
%hold on
%figure(2)
%plot(T2,Y2(:,1)/0.015,'b')
%hold on
%figure(3)
%plot(T2,Y2(:,3)/0.015,'G')
%hold on
figure(4)
plot(T2,Y2(:,1)/0.015,'-r',T2,Y2(:,3)/0.015,'--b')
xlabel('Time','FontName','Times New Roman','FontSize',20)
ylabel('w_1,w_2','FontName','Times New Roman','FontSize',20)
legend('uncontrol','control')
hold on

子程序：
function dy=control(t,y,U)
global kesai omega0 m rou D a1 a3 K C %%主系统基本参数
global epsilong_T s_33 C_pzt C_pzt1 K_pzt K_pzt1 k0 d_33 s_0 l_0 theta L R %压电材料参数
dy=[y(2);
    (rou*D*U*a1*y(2)*0.5+rou*D*a3*y(2)^3./U*0.5-C*y(2)-K*y(1))./m-K_pzt1*y(1)./m+theta*y(3)./m;
    y(4);
    -R*y(4)./L-y(3)./(L*C_pzt1)+theta*y(1)./L];
end
回复此楼

» 猜你喜欢

» 本主题相关商家推荐: (我也要在这里推广)

» 本主题相关价值贴推荐，对您同样有帮助:
I am what i am
1楼2015-11-09 15:56:56
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

getengqing

木虫 (正式写手)
这么一堆，谁会仔细看啊，建议楼主问详细一点的问题···
赞一下 回复此楼
一起交流学习/分享优秀资源
2楼2015-11-09 16:49:24
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

hustboshao

铜虫 (小有名气)
引用回帖:
2楼: Originally posted by getengqing at 2015-11-09 16:49:24
这么一堆，谁会仔细看啊，建议楼主问详细一点的问题···

不好意思，参数写的有点多，主程序其实就是一个ODE45函数，很简单的，但是不知道为什么算不出来
赞一下 回复此楼
I am what i am
3楼2015-11-09 16:59:43
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

lsjc77

铜虫 (初入文坛)

【答案】应助回帖

★ ★ ★ ★ ★ ★ ★ ★ ★ ★
感谢参与，应助指数 +1
hustboshao: 金币+10, ★★★很有帮助 2015-11-12 10:15:42
换成ode23s或者15s，没仔细看你的程序，但是发现很多很小的数值，感觉你的问题应该是刚性的，刚性方程不能用ode45解。
还有一个建议就是无量纲化，让你要解的参数在 0 1 之间，或者至少在一个数量级上。
不过你要解的方程很少，如果只是用ode45时间长，那换成ode23s/15s应该会很快给你一个结果。
赞一下 回复此楼

» 本帖已获得的红花（最新10朵）

hustboshao
4楼2015-11-10 14:54:40
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

hustboshao

铜虫 (小有名气)
送红花一朵
引用回帖:
4楼: Originally posted by lsjc77 at 2015-11-10 14:54:40
换成ode23s或者15s，没仔细看你的程序，但是发现很多很小的数值，感觉你的问题应该是刚性的，刚性方程不能用ode45解。
还有一个建议就是无量纲化，让你要解的参数在 0 1 之间，或者至少在一个数量级上。
不过你要 ...

无量纲化是个好办法，谢啦
赞一下 回复此楼
I am what i am
5楼2015-11-12 10:16:14
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 hustboshao 的主题更新
51/1返回列表
普通表情 龙 兔 虎 猫
信息提示
关闭
请填处理意见
关闭