Re: [討論] 積分上限

[ejialan (eji)]

※ 引述《hsuanboy (壓力)》之銘言：
: f(x)=0.000000000001602*x^5+0.000000003780060*x^4-0.000000034210106*x^3
:      -0.000166783274301*x^2+0.000174904487267*x+1.103243573568212
:   th1
:  ∫  f(x)dx
:   0                1
: ------------- = --------- ==> 可求出th1
:   90              100
:  ∫  f(x)dx
:   0
:   th2
:  ∫  f(x)dx
:   th1               1
: ------------- = --------- ==> 可求出th2
:   90               100
:  ∫  f(x)dx
:   0
:   th3
:  ∫  f(x)dx
:   th2               1
: ------------- = --------- ==> 可求出th3
:   90               100
:  ∫  f(x)dx
:   0
: 以上是我想表達的形式 第一次的積分下限初始值是0 每次求出積分上限值之後
: 會在帶入到下一行的積分下限去做運算 之後又會再求出積分上限值
: 依此類推 做100次
: 下面是我自己寫的程式 程式有跑出來 但是最後得到角度很怪
: 再做積分時 我是用徑度去做就算 算出上限的值都為勁度
: for迴圈的最後有再換成角度的形式 角度是以g表示
: 程式如下:
: syms x theatai
: th(1)=0;
: f=0.000000000001602*x^5+0.000000003780060*x^4-0.000000034210106*x^3-0.000166783274301*x^2+0.000174904487267*x+1.103243573568212;
: area=int(f,0,pi/2);
: integralf=@(s,t) int(f,x,s,t);
: for i=1:100
: sol = solve(integralf( th(i),theatai )-area*1/100);
: d = eval(sol);
: e = d(find( 0<d & d<pi/2 ));
: th(i+1) = e;
: g(i) = e*(180/pi);
: end

針對你的函數是多項式

matlab有幾個函數是針對多項式寫的

在定義多項式時只要把係數降羃排成列向量

之後可以用polyval計算多項式的值

polyint可以算反導函數的係數

計算會比使用符號運算快

另外solve也可改成fsolve或fzero

避免解出多解還要挑的情況

以下是照你的程式改寫給你參考

th=zeros(1,100);
a=[0.000000000001602
    0.000000003780060
    -0.000000034210106
    -0.000166783274301
    0.000174904487267
    1.103243573568212]';
ai=polyint(a);
area=polyval(ai,pi/2)-polyval(ai,0);
integralf=@(s,t) polyval(ai,t)-polyval(ai,s)-area/100;
for i=1:100
sol = fzero(@(t) integralf(th(i),t),th(i));
th(i+1) = sol;
g(i) = sol*(180/pi);
end

角度差距是先遞減再遞增

因為你的函數圖形是先遞增再遞減

可以plot(diff(g))來看

: 奇怪的問題點在 : 實際用手去計算每個角度之間的間距值應該是會越來越大才對
:                  但是我算完後 每個角度的間距都是差0.9度 程式哪邊有寫錯呢

--
※ 發信站: 批踢踢實業坊(ptt.cc) 
◆ From: 140.121.146.175
→ hsuanboy:有個問題想問一下 我用figure;plot(g)去畫圖來看
→ hsuanboy:跟用plot(diff(g))來看 不一樣的地放是在哪邊呢 另外
→ hsuanboy:這個函數本來應該是 -90~90度之間 中間y值最高
→ hsuanboy:但是我只想取0~90度的範圍就好 函數應該不是為遞減形式嗎
→ ejialan:diff是把離散數據做差分 所以diff(g)可以看出你的角度差
→ ejialan:所以你的f(x) 引數到底是徑度還是弧度?因為你的說明和
→ ejialan:程式寫的不一致
→ hsuanboy:真的很謝謝你喔
→ ejialan:上面打錯應該是徑度還是角度 如果是角度的話積分上下限
→ ejialan:就要改成用角度
→ hsuanboy:每次再帶入函數的時候 都是用勁度的
→ hsuanboy:最後再把所有的勁度換成角度
→ hsuanboy:不好意思 剛發現我想表示的有點錯誤 在分母的積分上限90
→ hsuanboy:應該改成pi/2為勁度才對 再帶入函數 求出的上限都是勁度
→ hsuanboy:最後才把所有的勁度換成角度
推 hsuanboy:不好意思 說明有點錯誤