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

河北科技大学硕士学位研究生 2012——2013学年第二学期 《Matlab语言及应用》结课论文 学 院： 信息科学与工程学院 专 业： 电路与系统 姓 名： 张利超 学 号： S2012014011 经典龙格库塔法及变步长龙格库塔法 1.经典龙格库塔法及变步长龙格库塔法matlab代码 a.经典龙格库塔法文件: Rungkutta4.m function R=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;for i+1:N??? k1 = feval(f,x(i),y(i));????? k2 = feval(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);end b.变步长龙格库塔法文件:change_step_RK.m function change_step_RK(fun); p21=2^p-1; while x(end)xN; [x1,y1]=RK_f(fun,x(n),y(n),h); [x2,y2]=RK_f(fun,x(n),y(n),h/2); if abs(y1-y2)/p21AbsTol; while abs(y1-y2)/p21AbsTol; x2=x1; y2=y1; h=2*h; [x1,y1]=RK_f(fun,x(n),y(n),h); end else while abs(y1-y2)/p21AbsTol; x1=x2; y1=y2; h=h/2; [x2,y2]=RK_f(fun,x(n),y(n),h/2); end end [xa,ya]=RK_f(fun,h,x(n),y(n)); x(n+1)=xa; y(n+1)=ya; n=n+1; end plot(x,y,k); function [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.利用两种方法求解初值问题 0x1,h=0.01 a.经典龙格库塔法 function varargout=saxplaxliu(varargin) clc,clear x0=0;xn=1;y0=0;h=0.01; [y,x]=lgkt4j(x0,xn,y0,h); n=length(x); fprintf( i x(i) y(i) \n); for i=1:n fprintf(%2d %4.6f %4.6\n,i,x(i),y(i)); end plot(x,y,k) function z=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; for i=1:n-1 K1=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*K3); y1(i+1)=y1(i)+h/6*(K1+2*K2+2*K3+K4); end y=y1; 运行结果截取前20个 i x(i) y(i) 1 0.000000 0.000000 2 0.010000 0.010050 3 0.020000 0.020200 4 0.030000 0.0304