河北科技学硕士学位研究生
2012——2013学年第二学期
Matlab语言应结课文
学 院：
信息科学工程学院
专 业：
电路系统
姓 名：
学 号：
典龙格库塔法变步长龙格库塔法
matlab代码
a典龙格库塔法文件
function RRungkutta4(fabNya)
h (ba)N
x zeros(1N+1)
y zeros(1N+1)
x ahb
y(1) ya
for i+1N
k1 feval(fx(i)y(i))
k2 feval(fx(i)+h2y(i)+(h2)*k1)
k3 feval(fx(i)+h2y(i)+(h2)*k2)
k4 feval(fx(i)+hy(i)+h*k3)
y(i+1) y(i)+(h6)*(k1+2*k2+2*k3+k4)
end
b变步长龙格库塔法文件change_step_RKm
function change_step_RK(fun)
p212^p1
while x(end) [x1y1]RK_f(funx(n)y(n)h)
[x2y2]RK_f(funx(n)y(n)h2)
if abs(y1y2)p21 while abs(y1y2)p21 x2x1
y2y1
h2*h
[x1y1]RK_f(funx(n)y(n)h)
end
else
while abs(y1y2)p21>AbsTol
x1x2
y1y2
hh2
[x2y2]RK_f(funx(n)y(n)h2)
end
end
[xaya]RK_f(funhx(n)y(n))
x(n+1)xa
y(n+1)ya
nn+1
end
plot(xy＇k＇)
function [xaya]RK_f(funhxy)
k1fun(xy)
k2fun(x+h2y+h*k12)
k3fun(x+h2y+h*k22)
k4fun(x+hy+h*k3)
xax+h
yay+h*(k1+k2*2+2*k3+k4)6
两种方法求解初值问题
a典龙格库塔法
function varargoutsaxplaxliu(varargin)
clcclear
x00xn1y00h001
[yx]lgkt4j(x0xny0h)
nlength(x)
fprintf(＇ i x(i) y(i) ＼n＇)
for i1n
fprintf(＇2d 46＼n＇ix(i)y(i))
end
plot(xy＇k＇)
function zf(xy)
zcos(x)+sin(y)
function [yx]lgkt4j(x0xny0h)
xx0hxn
nlength(x)
y1x
y1(1)y0
for i1n1
K1f(x(i)y1(i))
K2f(x(i)+h2y1(i)+h2*K1)
K3f(x(i)+h2y1(i)+h2*K2)
K4f(x(i)+hy1(i)+h*K3)
y1(i+1)y1(i)+h6*(K1+2*K2+2*K3+K4)
end
yy1
运行结果截取前20
i x(i) y(i)
8
b．变步长龙格库塔法
function change_step_RK1(fun)
funinline(＇cos(x)+sin(y)＇)
AbsTol1e4
h001
p4
x00
xN1
y00
xx0
yy0
n1
p212^p1
while x(end) [x1y1]RK_f(funx(n)y(n)h)
[x2y2]RK_f(funx(n)y(n)h2)
if abs(y1y2)p21 while abs(y1y2)p21 x2x1
y2y1
h2*h
[x1y1]RK_f(funx(n)y(n)h)
end
else
while abs(y1y2)p21>AbsTol
x1x2
y1y2
hh2
[x2y2]RK_f(funx(n)y(n)h2)
end
end
[xaya]RK_f(funhx(n)y(n))
x(n+1)xa
y(n+1)ya
nn+1
end
fprintf(＇ i x(i) y(i)＼n＇)
for i1n
fprintf(＇2d ＼n＇ix(i)y(i))
end
plot(xy＇r＇)
function [xaya]RK_f(funhxy)
k1fun(xy)
k2fun(x+h2y+h*k12)
k3fun(x+h2y+h*k22)
k4fun(x+hy+h*k3)
xax+h
yay+h*(k1+k2*2+2*k3+k4)6
运行结果
i x(i) y(i)
运行结果变步长龙格库塔法产生误差|y（i+1）y（i）|典龙格库塔法运算量较单步步长越截断误差越着步长缩定求解范围完成步数增加 步数增加引起计算量增导致舍入误差严重积累推算验证
典四阶龙格库塔公式
节点 出发先h步长求出似值
公式局部截断误差
然步长折半取步长跨两步求似值跨步截断误差
较（）式（）式步长折半误差约减少
易列事估计式
样通检查步长折半前两次计算结果偏差
判定选步长否合适 具体说区分两种情况处理：
1 定精度 果 反复步长折半进行计算直 止时取终 作结果
2 果反复步长加倍直止时步长折半次结果
表面选择步长步计算量增加总体考虑合算
文档香网(httpswwwxiangdangnet)户传
《香当网》用户分享的内容，不代表《香当网》观点或立场，请自行判断内容的真实性和可靠性！
该内容是文档的文本内容，更好的格式请下载文档
河北科技学硕士学位研究生
2012——2013学年第二学期
Matlab语言应结课文
学 院：
信息科学工程学院
专 业：
电路系统
姓 名：
学 号：
典龙格库塔法变步长龙格库塔法
matlab代码
a典龙格库塔法文件
function RRungkutta4(fabNya)
h (ba)N
x zeros(1N+1)
y zeros(1N+1)
x ahb
y(1) ya
for i+1N
k1 feval(fx(i)y(i))
k2 feval(fx(i)+h2y(i)+(h2)*k1)
k3 feval(fx(i)+h2y(i)+(h2)*k2)
k4 feval(fx(i)+hy(i)+h*k3)
y(i+1) y(i)+(h6)*(k1+2*k2+2*k3+k4)
end
b变步长龙格库塔法文件change_step_RKm
function change_step_RK(fun)
p212^p1
while x(end)
[x2y2]RK_f(funx(n)y(n)h2)
if abs(y1y2)p21
y2y1
h2*h
[x1y1]RK_f(funx(n)y(n)h)
end
else
while abs(y1y2)p21>AbsTol
x1x2
y1y2
hh2
[x2y2]RK_f(funx(n)y(n)h2)
end
end
[xaya]RK_f(funhx(n)y(n))
x(n+1)xa
y(n+1)ya
nn+1
end
plot(xy＇k＇)
function [xaya]RK_f(funhxy)
k1fun(xy)
k2fun(x+h2y+h*k12)
k3fun(x+h2y+h*k22)
k4fun(x+hy+h*k3)
xax+h
yay+h*(k1+k2*2+2*k3+k4)6
两种方法求解初值问题
a典龙格库塔法
function varargoutsaxplaxliu(varargin)
clcclear
x00xn1y00h001
[yx]lgkt4j(x0xny0h)
nlength(x)
fprintf(＇ i x(i) y(i) ＼n＇)
for i1n
fprintf(＇2d 46＼n＇ix(i)y(i))
end
plot(xy＇k＇)
function zf(xy)
zcos(x)+sin(y)
function [yx]lgkt4j(x0xny0h)
xx0hxn
nlength(x)
y1x
y1(1)y0
for i1n1
K1f(x(i)y1(i))
K2f(x(i)+h2y1(i)+h2*K1)
K3f(x(i)+h2y1(i)+h2*K2)
K4f(x(i)+hy1(i)+h*K3)
y1(i+1)y1(i)+h6*(K1+2*K2+2*K3+K4)
end
yy1
运行结果截取前20
i x(i) y(i)
8
b．变步长龙格库塔法
function change_step_RK1(fun)
funinline(＇cos(x)+sin(y)＇)
AbsTol1e4
h001
p4
x00
xN1
y00
xx0
yy0
n1
p212^p1
while x(end)
[x2y2]RK_f(funx(n)y(n)h2)
if abs(y1y2)p21
y2y1
h2*h
[x1y1]RK_f(funx(n)y(n)h)
end
else
while abs(y1y2)p21>AbsTol
x1x2
y1y2
hh2
[x2y2]RK_f(funx(n)y(n)h2)
end
end
[xaya]RK_f(funhx(n)y(n))
x(n+1)xa
y(n+1)ya
nn+1
end
fprintf(＇ i x(i) y(i)＼n＇)
for i1n
fprintf(＇2d ＼n＇ix(i)y(i))
end
plot(xy＇r＇)
function [xaya]RK_f(funhxy)
k1fun(xy)
k2fun(x+h2y+h*k12)
k3fun(x+h2y+h*k22)
k4fun(x+hy+h*k3)
xax+h
yay+h*(k1+k2*2+2*k3+k4)6
运行结果
i x(i) y(i)
运行结果变步长龙格库塔法产生误差|y（i+1）y（i）|典龙格库塔法运算量较单步步长越截断误差越着步长缩定求解范围完成步数增加 步数增加引起计算量增导致舍入误差严重积累推算验证
典四阶龙格库塔公式
节点 出发先h步长求出似值
公式局部截断误差
然步长折半取步长跨两步求似值跨步截断误差
较（）式（）式步长折半误差约减少
易列事估计式
样通检查步长折半前两次计算结果偏差
判定选步长否合适 具体说区分两种情况处理：
1 定精度 果 反复步长折半进行计算直 止时取终 作结果
2 果反复步长加倍直止时步长折半次结果
表面选择步长步计算量增加总体考虑合算
文档香网(httpswwwxiangdangnet)户传
《香当网》用户分享的内容，不代表《香当网》观点或立场，请自行判断内容的真实性和可靠性！
该内容是文档的文本内容，更好的格式请下载文档
