matlab 龍格庫(kù)塔法 變步長(zhǎng)龍格庫(kù)塔法

《matlab 龍格庫(kù)塔法 變步長(zhǎng)龍格庫(kù)塔法》由會(huì)員上傳分享，免費(fèi)在線閱讀，更多相關(guān)內(nèi)容在行業(yè)資料-天天文庫(kù)。

1、河北科技大學(xué)碩士學(xué)位研究生2012——2013學(xué)年第二學(xué)期《Matlab語言及應(yīng)用》結(jié)課論文學(xué)院：信息科學(xué)與工程學(xué)院專業(yè)：電路與系統(tǒng)姓名：張利超學(xué)號(hào)：S2012014011經(jīng)典龍格庫(kù)塔法及變步長(zhǎng)龍格庫(kù)塔法1.經(jīng)典龍格庫(kù)塔法及變步長(zhǎng)龍格庫(kù)塔法matlab代碼a.經(jīng)典龍格庫(kù)塔法文件:Rungkutta4.mfunctionR=Rungkutta4(f,a,b,N,ya)h=(b-a)/N;x=zeros(1,N+1);y=zeros(1,N+1);x=a:h:b;y(1)=ya;fori+1:N???k1=feval(f,x(i),y(i));?????k2=f

2、eval(f,x(i)+h/2,y(i)+(h/2)*k1);?????k3=feval(f,x(i)+h/2,y(i)+(h/2)*k2);?????k4=feval(f,x(i)+h,y(i)+h*k3);?????y(i+1)=y(i)+(h/6)*(k1+2*k2+2*k3+k4);endb.變步長(zhǎng)龍格庫(kù)塔法文件:change_step_RK.mfunctionchange_step_RK(fun);p21=2^p-1;whilex(end)

3、,y(n),h/2);ifabs(y1-y2)/p21AbsTol;x1=x2;y1=y2;h=h/2;[x2,y2]=RK_f(fun,x(n),y(n),h/2);endend[xa,ya]=RK_f(fun,h,x(n),y(n));x(n+1)=xa;y(n+1)=ya;n=n+1;endplot(x,y,'k');function[

4、xa,ya]=RK_f(fun,h,x,y);k1=fun(x,y);k2=fun(x+h/2,y+h*k1/2);k3=fun(x+h/2,y+h*k2/2);k4=fun(x+h,y+h*k3);xa=x+h;ya=y+h*(k1+k2*2+2*k3+k4)/6;2.利用兩種方法求解初值問題0

5、('ix(i)y(i)');fori=1:nfprintf('%2d%4.6f%4.6',i,x(i),y(i));endplot(x,y,'k')functionz=f(x,y)z=cos(x)+sin(y);function[y,x]=lgkt4j(x0,xn,y0,h)x=x0:h:xn;n=length(x);y1=x;y1(1)=y0;fori=1:n-1K1=f(x(i),y1(i));K2=f(x(i)+h/2,y1(i)+h/2*K1);K3=f(x(i)+h/2,y1(i)+h/2*K2);K4=f(x(i)+h,y1(i)+h*K

6、3);y1(i+1)=y1(i)+h/6*(K1+2*K2+2*K3+K4);endy=y1;運(yùn)行結(jié)果截取前20個(gè)ix(i)y(i)10.0000000.00000020.0100000.01005030.0200000.02020040.0300000.03045050.0400000.04080060.0500000.05125070.0600000.06179980.0700000.07244990.0800000.083198100.0900000.094047110.1000000.104995120.1100000.116043130.12000

7、00.127190140.1300000.138436150.1400000.149781160.1500000.161225170.1600000.172767180.1700000.184407190.1800000.196146200.1900000.207982b．變步長(zhǎng)龍格庫(kù)塔法functionchange_step_RK1(fun)fun=inline('cos(x)+sin(y)');AbsTol=1e-4;h=0.01;p=4;x0=0;xN=1;y0=0;x=x0;y=y0;n=1;p21=2^p-1;whilex(end)

8、y1]=RK_f(fun,x(n),y(n),h);