四阶RK(龙格库塔)方法求初值问题的数值解

用四阶RK方法求初值问题的数值解。取步长h=0.2, 计算结果取8位小数。
语言不限,要求上机调试通过。主要是我对这个函数实在不熟悉。。。
y'=1-2xy/(1+x*x) (0<=x<=2)
y(0)=0
1楼这个答案我原原本本的百度到过,但是不知道如何把上面这个问题做相应的修改,得到正确答案,所以才提问的。希望能给一个对应的解答,得到8位小数的解,而不是一个图,谢谢。

第1个回答  2010-12-08
f=inline('-y+x+1','x','y'); %微分方程的右边项
dx=0.05; %x方向步长
xleft=0; %区域的左边界
xright=10; %区域的右边界
xx=xleft:dx:xright; %一系列离散的点
n=length(xx); %点的个数

y0=1;

%%(1)欧拉法
Euler=y0;
for i=2:n
Euler(i)=Euler(i-1)+dx*f(xx(i-1),Euler(i-1));
end

%%(2)龙格库塔法
RK=y0;
for i=2:n
k1=f(xx(i-1),RK(i-1));
k2=f(xx(i-1)+dx/2,RK(i-1)+k1*dx/2);
k3=f(xx(i-1)+dx/2,RK(i-1)+k2*dx/2);
k4=f(xx(i-1)+dx,RK(i-1)+k3*dx);
RK(i)=RK(i-1)+dx*(k1+2*k2+2*k3+k4)/6;
end

%%Euler和Rk法结果比较
plot(xx,Euler,xx,RK)
hold on

%精确解用作图
syms x
rightsolve=dsolve('Dy=-y+x+1','y(0)=1','x');%求出解析解
rightdata=subs(rightsolve,xx);%将xx代入解析解,得到解析解对应的数值
plot(xx,rightdata,'r*')
legend('Euler','Runge-Kutta','analytic')