河北科技大学硕士学位研究生
2012——2013学年第二学期
《Matlab语言及应用》结课论文
学 院：
信息科学与工程学院
专 业：
电路与系统
姓 名：
张利超
学 号：
S
经典龙格库塔法及变步长龙格库塔法
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 0.0
发布者：全栈程序员-站长，转载请注明出处：https://javaforall.net/225342.html原文链接：https://javaforall.net
![按键精灵2014微信群内给好友轮流发私信[通俗易懂]](https://javaforall.net/wp-content/uploads/2020/11/2020110817443450-480x300.jpg)
