matlab编的4阶龙格库塔法解微分方程的程序
2010-03—10 20：16
function varargout=saxplaxliu(varargin)
clc，clear
x0=0；xn=1。2；y0=1；h=0。1;
[y,x］=lgkt4j（x0,xn，y0,h)；
n=length（x)；
fprintf（’ i  x（i)  y(i)\n’)；
for i=1:n
    fprintf('%2d %4.4f ％4.4f\n’,i，x(i），y（i))；
end
function z=f（x，y）
  z=-2*x*y^2；
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；
  printf输出格式matlab
结果：
i  x（i）  y（i)
1 0.0000 1。0000
2 0。1000 0.9901
3 0。2000 0.9615
4 0.3000 0。9174
5 0.4000 0。8621
6 0.5000 0。8000
7 0。6000 0。7353
8 0。7000 0.6711
9 0。8000 0。6098
10 0。9000 0。5525
11 1.0000 0.5000
12 1.1000 0.4525
13 1.2000 0.4098
龙格库塔法
一、基本原理：
龙格—库塔（Runge—Kutta)方法是一种在工程上应用广泛的高精度单步算法。由于此算法精度高,采取措施对误差进行抑制，所以其实现原理也较复杂。该算法是构建在数学支持的基础之上的。对于一阶精度的欧拉公式有：
yi+1=yi+h*K1
K1=f（xi,yi)
当用点xi处的斜率近似值K1与右端点xi+1处的斜率K2的算术平均值作为平均斜率K*的近似值，那么就会得到二阶精度的改进欧拉公式：
yi+1=yi+h＊（ K1+ K2)/2
K1=f（xi,yi)
K2=f(xi+h，yi+h*K1)
依次类推，如果在区间［xi，xi+1］内多预估几个点上的斜率值K1、K2、……Km,并用他们的加权平均数作为平均斜率K*的近似值，显然能构造出具有很高精度的高阶计算公式。经数学推导、求解，可以得出四阶龙格－库塔公式，也就是在工程中应用广泛的经典龙格－库塔算法：
yi+1=yi+h＊（ K1+ 2＊K2 +2＊K3+ K4)/6
K1=f（xi，yi)
K2=f（xi+h/2，yi+h＊K1/2）
K3=f(xi+h/2,yi+h＊K2/2）
K4=f(xi+h，yi+h*K3）
通常所说的龙格-库塔法是指四阶而言的，我们可以仿二阶、三阶的情形推导出常用的标准四阶龙格-库塔法公式
（1）
计算公式（1）的局部截断误差是 。
龙格—库塔法具有精度高,收敛，稳定(在一定条件下)，计算过程中可以改变步长,不需要计算高阶导数等优点，但仍需计算 在一些点上的值,如四阶龙格—库塔法每计算一步需要计算四次 的值，这给实际计算带来一定的复杂性,因此，多用来计算“表头”。
二、小程序
#include〈stdio.h>
＃include〈math.h>
#define f(x,y) （—1*（x)＊(y)＊（y))
void main(void)
{
double a,b,x0,y0,k1，k2,k3,k4,h；
int n,i；
printf("input a，b,x0，y0，n:”)；
scanf（”％lf％lf％lf%lf％d",&a，&b，&x0,&y0,&n)；
printf（”x0\ty0\tk1\tk2\tk3\tk4\n");
for（h=（b—a)/n，i=0;i！=n；i++)
｛
k1=f(x0,y0）;
k2=f（x0+h/2,y0+k1*h/2）;
k3=f（x0+h/2,y0+k2*h/2）;
k4=f(x0+h,y0+h*k3）;
printf(”％lf\t％lf\t"，x0，y0）;
printf（"％lf\t％lf\t",k1,k2)；
printf(”%lf\t％lf\n",k3,k4);
y0+=h*（k1+2*k2+2*k3+k4)/6；
x0+=h；
｝
printf（"xn=%lf\tyn=％lf\n”，x0,y0);
｝
运行结果：
input a，b，x0，y0,n：0 5 0 2 20
x0      y0      k1      k2      k3      k4
0。000000        2。000000        -0。000000      —0。500000      —0.469238
—0.886131
0.250000        1。882308        -0。885771      -1。176945      —1.129082
-1。280060
0.500000        1。599896        —1。279834      —1。295851      -1。292250
-1。222728
0.750000        1.279948        -1。228700      —1.110102      —1.139515
—0.990162
1。000000        1.000027        —1。000054      -0。861368      —0。895837
—0.752852
1.250000        0。780556        —0.761584      -0。645858      -0.673410
-0.562189
1.500000        0。615459        -0。568185      —0.481668      -0.500993
—0.420537
1.750000        0.492374        —0。424257      —0.361915      -0.374868
-0。317855
2.000000        0。400054        -0。320087      —0。275466      —0.284067
—0。243598
2.250000        0。329940        —0.244935      -0。212786      —0.218538
-0。189482
2。500000        0.275895        -0.190295      —0。166841      —0。170744
—0。149563
2.750000        0。233602        —0。150068      —0.132704      —0.135399
—0.119703
3.000000        0.200020        -0。120024      -0。106973      -0。108868
—0.097048
3。250000        0.172989        -0.097256      —0。087300      —0.088657
—0.079618
3.500000        0.150956        -0.079757      —0。072054      -0.073042
—0。066030
3.750000        0.132790        —0.066124      —0。060087      —0.060818
-0。055305
4.000000        0.117655        —0.055371      -0。050580      —0.051129
—0。046743
4.250000        0。104924        -0.046789      —0.042945      —0。043363
—0.039833
4。500000        0.094123        —0.039866      -0.036750      —0。037072
-0。034202
4.750000        0。084885        —0.034226      —0.031675      -0。031926
-0.029571
xn=5。000000     yn=0.076927

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。