变步长龙格库塔法matlab代码,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)

[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)/p21

while abs(y1-y2)/p21

x2=x1;

y2=y1;

h=2*h;

[x1,y1]=RK_f(fun,x(n),y(n),h);

end

else

while abs(y1-y2)/p21>AbsTol;

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.利用两种方法求解初值问题

0

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