作者ejialan (eji)
看板MATLAB
標題Re: [討論] 積分上限
時間Thu Jun 20 11:00:53 2013
※ 引述《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
1F：→ hsuanboy:有個問題想問一下 我用figure;plot(g)去畫圖來看 06/20 11:31
2F：→ hsuanboy:跟用plot(diff(g))來看 不一樣的地放是在哪邊呢 另外 06/20 11:32
3F：→ hsuanboy:這個函數本來應該是 -90~90度之間 中間y值最高 06/20 11:34
4F：→ hsuanboy:但是我只想取0~90度的範圍就好 函數應該不是為遞減形式嗎 06/20 11:36
5F：→ ejialan:diff是把離散數據做差分 所以diff(g)可以看出你的角度差 06/20 11:46
6F：→ ejialan:所以你的f(x) 引數到底是徑度還是弧度?因為你的說明和 06/20 11:54
7F：→ ejialan:程式寫的不一致 06/20 11:54
8F：→ hsuanboy:真的很謝謝你喔 06/20 11:54
9F：→ ejialan:上面打錯應該是徑度還是角度 如果是角度的話積分上下限 06/20 11:57
10F：→ ejialan:就要改成用角度 06/20 11:58
11F：→ hsuanboy:每次再帶入函數的時候 都是用勁度的 06/20 12:12
12F：→ hsuanboy:最後再把所有的勁度換成角度 06/20 12:13
13F：→ hsuanboy:不好意思 剛發現我想表示的有點錯誤 在分母的積分上限90 06/20 12:17
14F：→ hsuanboy:應該改成pi/2為勁度才對 再帶入函數 求出的上限都是勁度 06/20 12:18
15F：→ hsuanboy:最後才把所有的勁度換成角度 06/20 12:18
16F：推 hsuanboy:不好意思 說明有點錯誤 06/20 12:29