- 积分
- 16831
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2014-9-25
- 最后登录
- 1970-1-1
|
发表于 2019-11-19 19:28:49
|
显示全部楼层
回帖奖励 +20 金钱
inline和符号计算相当低效,你还把它放在二重循环里,就不要期待能给你算出结果。
把符号计算的步骤放在循环外面,只计算一次。用匿名函数的技巧将计算结果传给循环,这样速度快了n个数量级。
- clc,clear all,close all
- %执行必要的符号计算
- syms x p s1 k1 s2 k2
- fun=p-((23/89)*(1-exp(-((x-86.0175)./s1).^k1))+(66/89)*(1-exp(-((x-86.0175)./s2).^k2)));
- dfun=diff(fun,x); %导数
- %技巧:将符号计算结果转化成匿名函数
- eval(['fun=@(x,p,s1,k1,s2,k2) ',vectorize(fun),';']);
- eval(['dfun=@(x,p,s1,k1,s2,k2) ',vectorize(dfun),';']);
- %{
- 匿名函数求值的示例用法
- x=86;p=3;s1=3;k1=2;s2=5;k2=2;
- fun(x,p,s1,k1,s2,k2)
- dfun(x,p,s1,k1,s2,k2)
- %}
- %你原来的程序。略有改动
- quan=[0.5:0.001:0.999]; k=length(quan);
- for i=1:k
- for j=1:10000
- p=quan(i);s1=n(i,1);k1=n(i,2);s2=t(i,1);k2=t(i,2);
- [x_star,index,it] = newton(fun,dfun,100,p,s1,k1,s2,k2);
- y(j,i)=x_star;
- end
- end
- return
- %newton函数,有改动
- function [x_star,index,it] = newton(fun,dfun,x, p,s1,k1,s2,k2)
- it_max=100 ;
- ep=1e-5;
- index=0;k=1;
- while k<it_max
- x1=x;
- f1= fun(x1,p,s1,k1,s2,k2); %匿名函数求值,速度相当快
- f2=dfun(x1,p,s1,k1,s2,k2);
- x=x1-f1/f2;
- if abs(x-x1)<ep
- index=1;break;
- end
- k=k+1;
- end
- x_star=x;it=k;
- end
复制代码
|
|