24小时热门版块排行榜    

查看: 1673  |  回复: 3

赵小夭annie

铁虫 (初入文坛)

[求助] 怎样在fsolve中对初值进行循环并缩短运行时间?

本人在解四个非线性方程,每一个非线性方程需要给定对应初值才能解,本人给初值写了四个循环,分别对应的是t0i,t1i0,t2i0,t3i0。要求是当非线性方程解出的结果满足一定条件时则这组初值可选,问题在于四个for循环计算时间很长,请问大家有何方法可以改进?我用了fsolve解初值,把初值设置成了变量。以下为可运行程序,但是运行时间很长。

% Find the initial condition
clc
clear
tic
global t1i C1 C2 a d d1 d2 omega w yi t0i  t2i t3i t4i t1i0 t2i0 t3i0 t4i0 j d3 d4 C3 C4

a=20;
c=100;
d=0.5;
w=sqrt(c-d^2);
j=0;
m=0;
n=1;
A=zeros(10,n);
omega=4.8;
T=2*pi/omega;
t4i0=t0i+T;
d1=(c-omega^2)/((c-omega^2)^2+(2*d*omega)^2);
d2=(2*d*omega)/((c-omega^2)^2+(2*d*omega)^2);
d3=(-1)/(omega^2+4*d^2);
d4=(2*d)/((omega)*(omega^2+4*d^2));



% 以下为循环

for yi=5:0.1:13
    for t0i=0(2*pi/omega)/500)2*pi/omega)
            for t1i0=(1/5)*T+t0i:11/4)*T+t0i
                for t2i0=(2/5)*T+t0i:11/2)*T+t0i
                    for t3i0=(3/5)*T+t0i:14/5)*T+t0i
                        
                        
                           

C1=-a*(d1*cos(omega*t0i)+d2*sin(omega*t0i));
C2=(1/w)*(yi-a*((d2*d-d1*omega)*sin(omega*t0i)+(d1*d+d2*omega)*cos(omega*t0i)));

t1i=fsolve(@(t1i) (C1*cos(w*(t1i-t0i))+C2*sin(w*(t1i-t0i)))*exp(-d*(t1i-t0i))+a*(d1*cos(omega*t1i)+d2*sin(omega*t1i)),t1i0);
x2=((C2*w-C1*d)*cos(w*(t1i-t0i))-(C1*w+C2*d)*sin(w*(t1i-t0i)))*exp((-d)*(t1i-t0i))-a*omega*(d1*sin(omega*t1i)-d2*cos(omega*t1i));





C3=(1/(-2*d))*(x2+(a*omega)*(d3*sin(omega*t1i)-d4*cos(omega*t1i)));
C4=(1/(2*d))*(x2+2*d*1-(a/omega)*sin(omega*t1i));

t2i=fsolve(@(t2i) C3*exp(-2*d*(t2i-t1i))+C4+a*(d3*cos(omega*t2i)+d4*sin(omega*t2i))+1,t2i0);
x4=-2*d*C3*exp(-2*d*(t2i-t1i))-a*omega*(d3*sin(omega*t2i)-d4*cos(omega*t2i));





C1=-a*(d1*cos(omega*t2i)+d2*sin(omega*t2i));
C2=(1/w)*(x4-a*((d2*d-d1*omega)*sin(omega*t2i)+(d1*d+d2*omega)*cos(omega*t2i)));

t3i=fsolve(@(t3i) (C1*cos(w*(t3i-t2i))+C2*sin(w*(t3i-t2i)))*exp(-d*(t3i-t2i))+a*(d1*cos(omega*t3i)+d2*sin(omega*t3i)),t3i0);
x6=((C2*w-C1*d)*cos(w*(t3i-t2i))-(C1*w+C2*d)*sin(w*(t3i-t2i)))*exp((-d)*(t3i-t2i))-a*omega*(d1*sin(omega*t3i)-d2*cos(omega*t3i));
   




C3=(1/(-2*d))*(x6+(a*omega)*(d3*sin(omega*t3i)-d4*cos(omega*t3i)));
C4=(1/(2*d))*(x6-2*d-(a/omega)*sin(omega*t3i));

t4i=fsolve(@(t4i) C3*exp(-2*d*(t4i-t3i))+C4+a*(d3*cos(omega*t4i)+d4*sin(omega*t4i))-1,t0i+T);
x8=-2*d*C3*exp(-2*d*(t4i-t3i))-a*omega*(d3*sin(omega*t4i)-d4*cos(omega*t4i));


%选初值的条件

  if  abs(t4i-T-t0i)<0.1 && abs(x8-yi)<0.1 && t4i>t3i && t3i>t2i && t2i>t1i && x2<0 && x4<0 && x6>0 && x8>0
      figure(1)
      hold on
      axis([0 8 0 20])
      plot(omega,yi,'o')
      A(:,n)=[t0i;yi;t1i;x2;t2i;x4;t3i;x6;t4i;x8];
      n=n+1;
  end
  j=j+1;
  j
                    end
                end
            end
    end
end

toc
load chirp
sound(y,Fs)
回复此楼
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

竹一拿下

铜虫 (正式写手)

迭代啊。干啥要用for循环啊,这个最慢了

发自小木虫Android客户端
2楼2018-09-12 07:49:16
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

赵小夭annie

铁虫 (初入文坛)

引用回帖:
2楼: Originally posted by 竹一拿下 at 2018-09-12 07:49:16
迭代啊。干啥要用for循环啊,这个最慢了

意思是用牛顿迭代法解方程?那也是需要给循环的呀?
3楼2018-09-12 08:12:40
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

竹一拿下

铜虫 (正式写手)

4楼2018-09-12 08:55:39
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 赵小夭annie 的主题更新
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[硕博家园] 售SCI文章,我:8O.5.5.1O.54,科目全,可十急 +4 gy1nBQXYQJqL 2026-08-29 4/200 2026-08-30 01:49 by ZPa0EcMwuECS
[找工作] 售SCI文章,我:8O5.5.1.O.54,科目齐全,可+急 +4 gy1nBQXYQJqL 2026-08-29 6/300 2026-08-30 01:47 by ZPa0EcMwuECS
[找工作] 售SCI一区T0P文章,我:8O.55.1.O.54,科目全,可伽急 +3 gy1nBQXYQJqL 2026-08-29 6/300 2026-08-30 01:36 by ZPa0EcMwuECS
[公派出国] 售SCI文章,我:8O.5.5.1O.54,科目全,可十急 +5 ASdOkHsho7FD 2026-08-28 7/350 2026-08-30 01:12 by ZPa0EcMwuECS
[考博] 售SCI一区T0P文章,我:8.O.55.1.O54,科目全,可伽急 +6 ASdOkHsho7FD 2026-08-28 9/450 2026-08-30 01:04 by ZPa0EcMwuECS
[考研] 售一区SCI文章T0P,我:8O.551.O54,科目全,可十急 +5 ASdOkHsho7FD 2026-08-28 7/350 2026-08-30 01:03 by ZPa0EcMwuECS
[硕博家园] 售SCI一区文章,我:8.O.55.1.O.54,科目齐全,可伽急 +6 ASdOkHsho7FD 2026-08-28 10/500 2026-08-30 00:53 by ZPa0EcMwuECS
[硕博家园] 售SCI一区T0P文章,我:8.O55.1.O.54,科目全,可十急 +6 ASdOkHsho7FD 2026-08-28 12/600 2026-08-30 00:31 by ZPa0EcMwuECS
[基金申请] 国社科申报系统有变化 +3 kynobel 2026-08-23 3/150 2026-08-29 21:57 by 余韵清
[论文投稿] 售SCI一区文章,我:8.O.551.O.5.4,科目全,可伽急 +3 gy1nBQXYQJqL 2026-08-29 3/150 2026-08-29 18:26 by 4FFAWE8HcgUD
[教师之家] 导师吐槽:我怎么摊上了这么个极品研究生! +9 苏东坡二世 2026-08-23 9/450 2026-08-29 14:41 by hustersqt
[基金申请] 2026年叶企孙基金 +4 bud_bud 2026-08-27 7/350 2026-08-29 07:23 by foolishmani
[基金申请] 基金系统什么内容也没有 30+4 winsaint 2026-08-27 9/450 2026-08-28 11:06 by maolC
[基金申请] 哪位高人中了,把查询到的截图贴出来让我看看,让我长长见识 +5 yuleib84 2026-08-26 6/300 2026-08-28 00:02 by yudaoqian88
[基金申请] 为什么 国际(地区)合作与交流项目 没有放榜? 10+3 majunge000 2026-08-26 11/550 2026-08-27 08:42 by 北京莱茵编辑
[基金申请] 出来了 +9 trojank 2026-08-26 9/450 2026-08-26 14:25 by 宝贝虫子
[基金申请] 为什么国自然不能直接公布 +4 bjdxyxy 2026-08-26 4/200 2026-08-26 13:12 by qingmu1201
[基金申请] 今天务委会开完了,明天出结果吗 +19 angus9576 2026-08-25 23/1150 2026-08-26 10:03 by zp519
[基金申请] 国合现在查不到了吗? +10 chengyan1220 2026-08-24 20/1000 2026-08-26 08:57 by peasantsprig
[基金申请] 没有任何消息-是不是就凉了 +9 图啦图啦 2026-08-24 10/500 2026-08-25 11:59 by 南海小哥
信息提示
请填处理意见