chap7微分與數(shù)值積分_第1頁
chap7微分與數(shù)值積分_第2頁
chap7微分與數(shù)值積分_第3頁
chap7微分與數(shù)值積分_第4頁
chap7微分與數(shù)值積分_第5頁
已閱讀5頁,還剩105頁未讀, 繼續(xù)免費閱讀

下載本文檔

版權(quán)說明:本文檔由用戶提供并上傳,收益歸屬內(nèi)容提供方,若內(nèi)容存在侵權(quán),請進(jìn)行舉報或認(rèn)領(lǐng)

文檔簡介

第七 數(shù)值微分與數(shù)值積7‐1:已測得某地大氣壓強隨高度變化的一組數(shù)據(jù)高度0壓強kgf1000米處的大氣壓強變化率。注釋:已知函數(shù)點的值,求在某些點的導(dǎo)數(shù)bI[f]b

f(x)dxF(b)Ffxsinx,fxex2x 我們是否可以找簡單方法計算下列積分11 11 x , xe2xdx dx 01

x22 x(2x)2

dx??x12345y468‐yfx未知只知道部分對應(yīng)點x12345y468注釋:通過函數(shù)點的值,求某個區(qū)間上的積分積分=已知數(shù)據(jù)點上函數(shù)值的線性§7.1x""x""y""BfABfAC

x各種一階差商f(x)f(xh)fhf(x)f(x)f(xhf(x)f(xh)f

向前差商公式D向后差商公式D中心差商公式Dff(xh)f(x)hf(x)23f(x)f(x)f(xh)f(x)hf(x)23f(x)f(x)向前:fxfxhfx)hf(

向后:fx

f(x)f(xh)

f(

f(xh)f(xh)

f( 中心:f 二階差商f''(x)f(xh)2f(x)f (x)f(xh)2f(x)f(xh)2f(4)注意例設(shè)f(x=ex分別取步長h102103,104106用向前差商公式和中心差商公式計算f(1.15)hxy向前差商中心差商上表結(jié)果說明當(dāng)步長由102縮小到104時,誤差在縮小.進(jìn)而縮小到106時,由于有效數(shù)字丟失,誤差反而增大。當(dāng)步長不太大時,中心差商公式精度優(yōu)于向前差商公式。以插值多項式的導(dǎo)數(shù)作為函數(shù)導(dǎo)數(shù)的近似,f(xi)n(xiBfABfAC

x插值函數(shù)的誤差公式f(n1)R(x)f(x)n(x) (n1)!插值型微分的誤差公R(x)f(x)(x)

(x)

f(n1)()

(n1)

對任意x[a,] 因(依賴于x)未知,故上式很難估計誤差,但若只求某個節(jié)點上的導(dǎo)數(shù)值,誤差可估計.因此,插值型求導(dǎo)公式通常用于求節(jié)點處導(dǎo)數(shù)的近似值。

R(

)f(

)n(xi)

(n

(xixjji1等距節(jié)點常用公式1.兩點公1

x0,x1x0f(x)

f(

)(xx1

f(x

(xx0 (

x1

(

x0f(x)L(x)

f(x1)f(x0

i R(x)f(x)L(x)f(0)(xx)hf(

R(x)f(x)L(x)f(1)(xx)hf(

ff(x)f(x1)f(x0)0hR(x)hf(1020f(x)f(x1)f(x0)1hR(x)hf( 21,(a, 2.三點公

xix0

(i二次插值公f(x)f(x

(xx1)(x

f(x

(xx0)(x

f(x

(xx0)(x0

x)(xx

1

x)(xx

2

x)(xx

f(x)(xx1)(xx2)f(x)(xx0)(xx2)f(x)(xx0)(x

插值公式求L(x)f(x)2x

f(x

2x

f(x)2x

因此f(x)L(x)f(x)2x0x1

f(x)2x0x0

f(x)2x0x0

f(x0

h2f

f(x2)

3f(x0)4f(x1)

f(x2f(x)L(x)f(x)2x1x1

f(x)2x1

f(x)2x1x0 hf(x)

f(x)

f(x2)f(x0

f(x)L(x)

f(x)2x2x1

f(x)2x2x0

f(x)2x2x0

f(x)

f(x)

f(x)

f(x0)4f(x1)3f

f(x f(xf(x)

3f(xi)4f(xi1)f(xi2 i f(xi)

f(xi2)4f(xi1)3f(xix i插值公式二次求導(dǎo)L(x)f(x f(x)2f(x

f(x)f(x0)2f(x1)f(x2 i0,1, xix0

(if(xf(x)13f(x)4f(x)f(x0012f(x)1f(x)f(x102f(x)1f(x)4f(x)3f(x2012(x0) (002R(x)2

f(3)( 2R(x)

f(3)( i(a,

(i0,1,ff(x)(x)1f(x)2f(x)f(xi2i012設(shè)Sx)fx)的三次樣條插值函數(shù),由三次樣條f(x)S(且有誤差估

xR(k)(x)f(k)(x)S(k)(x)O(h4k (kf(x)f(xh)f(f(x)hf(x)f(xh

f(x)f(x)13f(x)4f(x)f(x0012f(x)1f(x)f(x102f(x)1f(x)4f(x)3f(x2012f''(x)f''(x)f(xh)2f(x)f(x§7.2Newton-Cotesfx)的原函數(shù)不存在或計算繁瑣,或fx)只有b離散數(shù)據(jù)點,求I(f) f(x)dx的近似數(shù)值fn(x)bI(fn) fn(x)dx作為If)bEn(f) f(x)fn(x) n nbI(f)

f(x)dx

Ln(x)dxlk(x)dxf(xkk 插值型求積公式b

I(f)

f(x)dxAkf(xkkkAkk

lx)dx(k01n)

f無關(guān),稱為求積系求積公式的截斷誤

f(n1)(R(f)R(f)f(x)dxL(x)dx (

(n

xa (k0,1,",n),hbk b

x Ak

lk(x)dx

j

xxjxxa

j t

(1)nk(b

kk

hdt

k!(nk

0(tjkk

(1)nkjnk!(nk)!nn

(tnn

jk(ba)Ckb

n(I(f)

f(x)dx(ba)k

f(xkC(n)

(1)nk k!(nk)!n

(tjjbn階N-C公式的截斷bRn(f)

f(n1)((n

(xxj)dxnn

nh(n1)!nh

f(n1)()(tnnn12Cn12C(nk12123879038790n 梯形公I(xiàn)(f)baI(f)baf(a)4f(ab)f62I(f)f(x)dxbaf(a)fba2梯形公式

f(x)dxhf(x)f(x)b 2b

Simpson公式

f(x)dxhf(x)

4f(x)f(xb 3b

Simpsons3公式8bf(x)dx3hf(x)3f(x)3f(x)f(xb 8

衡量插值型求積公式的精度,用多項式次數(shù)作標(biāo)(一)求定義7.1

f(x)dxAkf(xkbnkbnbn f(x)dxAkf(xkbnk存在某個m1次多項式Pm1x),bb

AkPm1(xknkn則稱該求積公式具有m次代數(shù)精確度結(jié)論

f(x)dxAkf(xkbnkbn

具有m次代數(shù)精度n bxldxAxn

(l0,1,"ma k m

m

dx

k

Akxk由n次插值導(dǎo)出的求積公式至少具有n次代數(shù)精度bRn(f)b

f(x)dx

nkn

Akf(xk)

b(n n1(bnR(xl)n

(l0,1,",結(jié)論2n階NC公式至少具有2n1次代數(shù)精度 設(shè)2n1P2n1x

x2n1

"ax

P(2n1)

bb bxanh

nn

(

h2n

-

njn

(t+n-j)

h2n

-

t(t

1)(t222)"(t

n2)對稱區(qū)間上奇函數(shù)積分為 證bf(x)dxb

ba

f(a)4f(ab)fa

ba a 記I(f

f(

I(f)

f(a)4f )f I(x)1(b2

I~

ba(141)ba6ba(a2a2bb)1(b2I I(x2)1(b3a3)

ba

ab

2 I(x)

a

)b

(ba)I(x3)1(b4

baa34(ab)3b3

4a4)4 I(x4

I(x4)1(b5

4ba

44(ab)4b4I(x5

ba(5a44a3b6a2b24ab35b4) 例:確定系數(shù)使下列求積公式代數(shù)精確度盡可11

f(x)dx

f(1)

f(0)

ffx1xx2時公式等號成立 A1A0A1

A1

A0 A

A f(x)dx

f(1)

f(0)

f

11

fx1xx2精確成立fx)x3時,

A0 A當(dāng)f(x)x4 左邊2右邊2

定理7.2fxC(2abbb

f(x)dx

f1b證明:R(f1b

f(x)(xa)(x2bf(b

(xa)(xf()(ba)3 7.3fxC(4ab則SimpsonbR(f) f(x)dxb

baf(a)4f(ab)f(b)

(b 證明

f(4)() f(4)( (a, Simpson公式的代數(shù)精度為3尋找fx)插值多項式使得在a,abb點的值等于函數(shù)2值且是三次多項式。取fx)的三次Hermite插值多項式H3x)H3(a)f H(ab)f(a H(ab)f(a baf(a)4f(ab)f baH(a)4H(ab)H(b) bH( f(4)( af(x)H3(x) (xa)(x

)2(x bR(f) f(x)dxb

baf(a)4f(ab)f(b)

f(4)(

a f(x)dx

H3(

(xa)(x

)2(x2f(4)(

(xa)(x

a

)2(xb bxabtb f(4)( ba

)1(t1)t(tf(4)((b)5f(4)((b)522253(b

)

f(4)(

f(4 (a f(n1) f(n1)(R(f)ban1(x(n為奇數(shù)n(n(n2(n2)! ba(xa2(x(n為偶數(shù)其 n1(x)(xx0)(xx1)"(xxn (a, xka(k用 1",xn

n(k0,1,",

xn1II(f)baf(x)dx(ba)bkf(xk)En(fknb bl(

(k0,1,", ba 次數(shù);但是高次插值有Runge現(xiàn)象,不一定能提復(fù)復(fù)化求積基本思用分段低次多項式近似被積函數(shù)求積分近似復(fù)化梯形fx) a

(k0,1,",

yhbyhbnnbf(x)dxb

xk

f(

hf(x)f( k xkk

k0

kf(a)f(b)2

f(xk k bbf(x)dxhf(a)a2f(b)f(x)kkn如果fxC(2)a其截斷誤差RR(f) (ba)f(2T(a, h

RT(f)

f(x)dx

f(a)f(b)n2 kn2

f(xkh h

f

)]

h2(ba)

f(

(x, k

k

kh2

(ba)f(

(a,fhbn分段Simpson公hbny hx1

h

n

hba

a

(k0,1,",x2x2k

f(x)dx2h[f( 2k

)4f(

2k1)f(x2k44444 44444mbf(x)dxm

f(x)dx2hf

)4f(x )f(x

m hf(a)f(b) m 3 f(x2k f(x2k1 k k baf(x)dxh3mf(a)f(b)kf()kmf(2k)RR(f)sbh (a, h

Rs(f)

f(x)dx

3f(a)f(b)

f(

)4f(x2k1m m

5

kkf(4)(k

k f(4)()

b

k

b

kk

h k(x2k2,x2k

(a,

0

x21問積分區(qū)間要等分為多少份才能保證計算結(jié)14位有效數(shù)字解:0.3e

x0,1

1ex2 0 11042

f(x)ex

f(x)2xef(x)(4x22)ex22f(4)(x)(16x448x212)e

f(x)(8x312x)e

Rs(f

f(4)(

1800

f(4)(

104n421049

n

n8即可1ex2dx0

1 (14e1 38

e1)f(x)(4x22)e

f(x)(8x312x)ex2

fRT(f)2

f(

120 h21104 n n7.3.3逐次分半算法復(fù)化梯形公式誤差h h (f)

f

[f()

k

k

kbb

f(x)dx

h[22

(b)

注意到區(qū)間再次對1h RT(f) [f(b)f(a)]4RT(f2 122 I13I13TnI kks (f)sn

m590km5

f(4)(

h4

k

(4)(kbh4kb180

f(4)(x)dx

(f(b)f(a))注意到區(qū)間再次對f(b)f(a)h (f)2

(fSnS I1(I1(SnI 通常采取將區(qū)間不斷對分的方法即取n取n2m(m0,1,2,"), ba

xa (k0,1,",2m 2m m m

f(a)f(b)2

f(x

(m0,1,2," 2

k 12mf(a(2k1)hm m(m1,2,"k即Tm是前次計算的近似值Tm1的一半與新分點處函數(shù)2 2m之和乘以新步長h之和m

Tm稱為梯形值序列IR(fIR(f )13)Sm2 Sm2 3f(a)f(b)f(a2kh)mf(a(2khm(m1,2,"2

fxC(4ab2IS R(f,2

)1(

Sm122S S 例:計算橢 4

的周長,使結(jié)果具有5位有效數(shù)解x2cosy lds4 4sin2cos2d4

3sin2L

則 0

2

|e(I) 813sin213sin2T(12) T1Tf

)

1|TT|

T1T f(3)

8f(8 1|TT| T1T(f()f(3)f(5)f(7)) 13|T4T8| 1

184T842.42211l4T820203sin2 (f(0)4f()

f( ))2. (f(0)4f(

)2f()4f(3)f())

1|SS|

[f(0)4f

)2f )4f( 2f(4

)4f(5)2f(3)4f(7)

f(

1|SS

1104 l4I42.42211梯形公式算法輸入ab,f取m1,hba Th[f(a)f 取F0,對k1,22m1,Ff(a(2k1)hT1T 5若|TT|輸出I停止m1h TT0轉(zhuǎn)Simpson公式算法輸入ab,fx取Ff(a)f(b),Ff(ab), ba(F4F m hb4

3F303

k1,2,"

m1,

f(a(2k1)h)Sh(F2F4F |SS|15ISm1h F2F3 S 轉(zhuǎn) 計算圓周率的一種算法n nsinn 7∵sinxx1x31x51x 7

(1

(1

(1)6 7 逐次分半時粗細(xì)近似值的組合 1 2 2 (22

1(1615215

(1)n2

1(42nn3(232

1(1615215

(1)nnn(1)1 2n (2)1 2n 6II(h)a

a

"ahPk 其中01P2" II(h)ah2ah4ah6 問題O(hp1)提高到O(hp2)?

II(h)ah2ah4ah6 II(h/2) 4a216a364

則(2)4 3I[4I(h/2)I(h)]3a24

15a3h6I4I(h/2)I(h)c

ch6 I(h

(h

ch6" 22m (h/2)

2(m1)似Im(h) 22m 則IIm(h) 1對 II(h)a1

ahP2"a 1II(h)aP1 1

aP2

"aPk P1IP1I2k3k2a2

P1)hP2a

P1

"a

P1)hPk(I(h)P1I(h))I (1P1)

a

"a

"(1P1

kI

f( 記 T(k)I(bab

(k0,1,2," 由Euler-MaclaurinIT(k)a(bIT(k)a(ba)2a(ba)4"a(ba)2i012iT(k) 4

則 R(f,T(k))O((ba)4 T(kT(k) m m4mT(k1)T(km4mRf

(k))mm

)2m2m上式Tkm

(m12稱為Romberg求積公式m以T (m0,1,2,")作為積分近似值,稱為Romberg方法m m停止準(zhǔn)則:T(0)T( m 逐行自左向右計分別將求積區(qū)間2k等分(逐次對分 T(k

42T(k1)

43T(k1) T01T010 T

T(

42

43 T24T240

T51T5

T632T637 70

818

929

(k)是Simpson序列 2kmmmRomberg方法產(chǎn)生的序列(0當(dāng)m時收斂于Ifm且收斂速度快于復(fù)化Simpson方 一般公 I(f)b

f(x)dxAkf(xkkk求積系A(chǔ)kk

lx)dx(k01n)與f等距節(jié)點:N-C公 逐次分半梯形,Simpson 啟示:在節(jié)點數(shù)n固定時適當(dāng)選取xk}和求積系數(shù)Ak可以提高數(shù)值積分方法的近似程度 基本思想節(jié)點數(shù)n固定適當(dāng)選取xk}和求積系數(shù)Ak使求積公式

(x)f(x)dxbnkbn

f(xk具有最高的代數(shù)精確度其中x0為權(quán)函數(shù),求積系數(shù)Ak(k01n)與f無關(guān)。 n個節(jié)點的求積公式最高代數(shù)精確度是多大怎樣選取節(jié)點與求積系數(shù)使求積公式具有最高的代數(shù)精確度?bn(一)Gaussbn

a(x)f(x)dx

k

精度,則節(jié)點xk}和求積系數(shù)Ak}應(yīng)滿足方程A1A2"AnAxA "A A

A

"Axm b b

(x)xldx

(l0,1,…,此2n個未知數(shù)的m+1個方程可證,當(dāng)m12n時有解設(shè) (x)(x

x)2(xx)2"(x

x

b a(x)P2n(x)dxbn AkP2n(xk)k所以n2n定義7.2如果一組節(jié)點x1x2,xnab能使求積公 a(x)f(x)dxAkf(xkk具有2n1次代數(shù)精度則稱這組節(jié)點為Gauss點上述求積公式為帶權(quán)函數(shù)(x)的Gauss型求積公式II1f(x)dxf(1)f(133I1f(x)dx5f935)89f(0)59f35)兩 公式三 公式例11f(x)dxaf 0.6)bf(0)cf 11中的待定系數(shù)ab和c使其代數(shù)精度盡量高,并,1解:If

f(x)dx,I(f)af 0.6)bf(0)cf 取fx1x,x2使I(fI(f1I(x)

1dx I(1)ab111xdx I(x) 0.6a11

I(x2)

x2dx2 I(x2)0.6a0.6c要使公式具有盡可能高的代數(shù)精度ababc

ac9 a

c 2

b0.6a0.6c ac5b8時求積公式的代數(shù)精度最高 1f(x)dx5f90.6)8f(0)5f990.6I(x3)

x3dx

5

0.6)3

I(x I(x4)

x4dx5

~(x4)

59

0.62)0.40.60.6 I(x)1x

I(x) 0.6)5 0.6

(6II(x6)1x6dx ~x)5(0.630.63)2.16(6I 4例: 求積公式計算4

xe2xdx解:通過坐標(biāo)變換化為[-1,11令x2t14I4

xe2xdx

(4t4)e4t

f(t(4t4)e4t4兩 公式

4

4I1f(t)dtf

)f

)(4

3(4 3I1f(x)I1f(x)dxf(1)f(133

三 公式

4xe2xdx0I

f(t)dt5f 0.6)8f(0)5f 115(49

0.6)e4

8(4)e45(4

II1f(x)dx5f953)8f(0)5f9935)5 )8 )5 注釋:Gauss型求積公式由代數(shù)精確度描述 最高階,達(dá)到2n-求積節(jié)點和系數(shù)有特殊的選取 或者(二)Gauss 確定x1,x2,A1,A2,使得求積公1

f(x)dxA1f(x1)

f(x2具有最高代數(shù)精確度。解:Gauss求積公式可以達(dá)到3次代數(shù)精確度.取fx)為任意三次多項式,利用多項式除法,fx)可表成如下形式f(x)(a0a1x)(xx1)(xx2)(b0

f(x)dx

a1x)(xx1)(xx2)dx

A1f(x1)A2f(x21(011)如果求積公式具有3次代數(shù)精確度,則11(b0b1x)dxA1(b0b1x1)A2(b0b1x21 1(a0a1x)(xx1)(xx2)dx對任意實數(shù)a0,a1都成立1 (xx)(xx)dx1 x(xx1)(xx2)dx即 22xx 22

(x1x2)3由此解出x= 13 再取兩個特殊的f,求得A1和A2:比如取f11f21 dx2A1

331- 33

xdx0 其解為A1A21。故所求求積公式為1f(x)dxf1f13333

例如f(x)xx2xx 則f(x)dx

1x-x2x-

2dx11而

f(x1)f(x2)b求積公式b

(x)f(x)dx

nkn

fxk)nn(x)(xxknk與任意至多n1次多項式qx)在[ab]上關(guān)于權(quán)函數(shù)x)正交即b(x)q(

(x)dx 若xk}為Gaussba(x)f(x)dxb

nkn

Akf(xk為帶權(quán)函數(shù)x)的Gauss型求積公式n則對任意至多2n1Px)nba(x)P(x)dxb

k

AkP(xk若qx)n1次多項式bn則qx)n(x)至多2n1次多項式bn

(x)q(x)n(x)dx

k

Akq(

)n(

)

(x)(xxkkf(xxl(l01n1)A1A2"AnAxAx"A

(x)xldxbab

lnAnA

A

A

"

Axn1

b

b(x)xldx

n

(k1,2,",ka(x)f(x)dxk

k

Akf(

代數(shù)精度至少為n設(shè)fx)為任意次數(shù)不超過2n1則存在次數(shù)不超過的n1多項式qx),r(x)使f(x)q(x)n(x)r(bb(x)f(x)dxb

(x)q(

(x)dx

b(x)r(a Akr(xk)Akf(xk bkb

k即求積公式

(x)f(x)dxk

Akfxk)對次數(shù)不超2n1次的多項式精確成 求積公式

(x)f(x)dx

k

f(

)為Gauss公式,{xk}是Gauss點 證回定理6.3設(shè)nx)(n0,12是最高次項系數(shù)不為零的n條件:對任意n1次多項式qx都有bb

(x)q(

(x)dx (n1,2,"求求Gauss knk[ab]上關(guān)于權(quán)函數(shù)x)的正交項式系中的n次多項式的n個實根A b(x)l( A b(

n(

(xx)(x因為lkx)

nnj

xxxkx

是n1次多項式,所jk b

n( a(x)lk(x)dx

(

(xx)(x)dx

Allk(xl)注:

bk(x)l2(bk

lk l2(k

xxj

jj

xkxjnn b(x)l2(x)dx a

l

Al2(x) 結(jié)論f(x)C(2n)a,b,則(a, 使bR(f)b

(x)f(x)dx

nkn

Akf(xkf(2n)((2n)!

b(x)2(x 證明:構(gòu)造過節(jié)點x1,x2xn上f(x)的2n-1階f2n( (x) (a, R(f)

(x)f(x)dxk

Akf(xk

a(x)f(x)dx

k

AkH(xk b(x)f(x)dxb(x)H(

f(2n)(x)bb

(x)2(nf(2n)()bn

n(n

2(

(a, fxC(ab),Gaussn時bb

(x)l(x)dxbkabk

kbakb

nkknk

(nnkn

f(xk)

nkn

Ak

(x

nkkknkk

f(x)f(x k

Akk1kkk

f(xk)

(xbb

(x)dx1k

f(xk)

(x

[1,1上帶權(quán)函數(shù)x1Gauss11

f(x)dxAkf(xknnGauss點xk}n次Legendre多項nPx) x21)nn個零點n 2nn!dxn 求積系

1,2,",n) (1x2)P

(2n

b求積分b

f(

做變

xbatba 11b f(x)dxb

b2

f(batba kbak

Af(ba

ba k

II1f(x)dxf(1)f(133I f(x)dx5f135)89f(0)59f935)11的Gauss型求積公1f(x1xndxAkf(xkkGauss點 cos2k1(k1,2,",n)為n次Chebyshev多項kTnxcos(narccosx)n個零

(k1,2,", f(

2k1x1x

dx

nkn

f(cos 插值余:R(f) f插值余

[0)上帶權(quán)函數(shù)x)ex的Gauss型求積公 f(x)dxAkf(

溫馨提示

  • 1. 本站所有資源如無特殊說明,都需要本地電腦安裝OFFICE2007和PDF閱讀器。圖紙軟件為CAD,CAXA,PROE,UG,SolidWorks等.壓縮文件請下載最新的WinRAR軟件解壓。
  • 2. 本站的文檔不包含任何第三方提供的附件圖紙等,如果需要附件,請聯(lián)系上傳者。文件的所有權(quán)益歸上傳用戶所有。
  • 3. 本站RAR壓縮包中若帶圖紙,網(wǎng)頁內(nèi)容里面會有圖紙預(yù)覽,若沒有圖紙預(yù)覽就沒有圖紙。
  • 4. 未經(jīng)權(quán)益所有人同意不得將文件中的內(nèi)容挪作商業(yè)或盈利用途。
  • 5. 人人文庫網(wǎng)僅提供信息存儲空間,僅對用戶上傳內(nèi)容的表現(xiàn)方式做保護(hù)處理,對用戶上傳分享的文檔內(nèi)容本身不做任何修改或編輯,并不能對任何下載內(nèi)容負(fù)責(zé)。
  • 6. 下載文件中如有侵權(quán)或不適當(dāng)內(nèi)容,請與我們聯(lián)系,我們立即糾正。
  • 7. 本站不保證下載資源的準(zhǔn)確性、安全性和完整性, 同時也不承擔(dān)用戶因使用這些下載資源對自己和他人造成任何形式的傷害或損失。

最新文檔

評論

0/150

提交評論