| 查看: 5108 | 回复: 4 | ||
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 |
回复此楼» 猜你喜欢
基金申请
已经有44人回复
-
CSC与新西兰维多利亚大学PhD奖学金项目
已经有0人回复
-
物理学I论文润色/翻译怎么收费?
已经有177人回复
-
新西兰Robinson研究所 招聘CSC公派访问人员
已经有0人回复
-
帮我的英语口语老师找学生
已经有0人回复
-
什么时候开奖？
已经有13人回复
-
散金币祈福
已经有94人回复
-
青基已中
已经有1人回复
-
散金币祈福
已经有87人回复
-
求助，如何提取ELK的rt-TDDFT在某一时刻的自旋密度分布和电子密度分布
已经有1人回复
» 本主题相关商家推荐: (我也要在这里推广)
» 本主题相关价值贴推荐，对您同样有帮助:
- 使用Multiwfn绘制原子轨道图形、研究原子壳层结构及相对论效应的影响
已经有17人回复
- 大家看看，我该如何投稿到哪个期刊？《从金庸的武学视角谈计算机编程思维的培养》
已经有10人回复
getengqing
木虫 (正式写手)
- 应助: 104 (高中生)
- 金币: 3759.8
- 散金: 30
- 红花: 5
- 帖子: 592
- 在线: 110.5小时
- 虫号: 2986702
- 注册: 2014-02-23
- 性别: GG
- 专业: 自然语言理解与机器翻译
2楼2015-11-09 16:49:24
hustboshao
铜虫 (小有名气)
- 应助: 0 (幼儿园)
- 金币: 1308.9
- 红花: 1
- 帖子: 137
- 在线: 78小时
- 虫号: 1567118
- 注册: 2012-01-06
- 性别: GG
- 专业: 固体力学
3楼2015-11-09 16:59:43
【答案】应助回帖
★ ★ ★ ★ ★ ★ ★ ★ ★ ★
感谢参与，应助指数 +1
hustboshao: 金币+10, ★★★很有帮助 2015-11-12 10:15:42
感谢参与，应助指数 +1
hustboshao: 金币+10, ★★★很有帮助 2015-11-12 10:15:42
|
换成ode23s或者15s，没仔细看你的程序，但是发现很多很小的数值，感觉你的问题应该是刚性的，刚性方程不能用ode45解。 还有一个建议就是无量纲化，让你要解的参数在 0 1 之间，或者至少在一个数量级上。 不过你要解的方程很少，如果只是用ode45时间长，那换成ode23s/15s应该会很快给你一个结果。 |
» 本帖已获得的红花（最新10朵）
hustboshao4楼2015-11-10 14:54:40
hustboshao
铜虫 (小有名气)
- 应助: 0 (幼儿园)
- 金币: 1308.9
- 红花: 1
- 帖子: 137
- 在线: 78小时
- 虫号: 1567118
- 注册: 2012-01-06
- 性别: GG
- 专业: 固体力学
5楼2015-11-12 10:16:14
30