matlab 龙格库塔法 变步长龙格库塔法

《matlab 龙格库塔法 变步长龙格库塔法》由会员上传分享，免费在线阅读，更多相关内容在行业资料-天天文库。

1、河北科技大学硕士学位研究生2012——2013学年第二学期《Matlab语言及应用》结课论文学院：信息科学与工程学院专业：电路与系统姓名：张利超学号：S2012014011经典龙格库塔法及变步长龙格库塔法1.经典龙格库塔法及变步长龙格库塔法matlab代码a.经典龙格库塔法文件: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.变步长龙格库塔法文件: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;运行结果截取前20个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．变步长龙格库塔法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);