利用MATLAB完成程序编写

编写二阶龙格-库塔方法的程序,并求解西门的初值问题,写出程序源代码,输出计算结果.

大概就是这样的了，时间紧急，谢谢了！！！

function z=lkfun(x,y) %f(x,y)部分，可以根据具体的函数修改

z=-0.9*y/(1+2*x);

%求解函数

%四阶龙格库塔常微分方程数值解MATLAB编程

function [y,x]=LK(a,b,y0,N)

%a,b表示数值解的区间

%y0表示初始值

%N表示解的空间密度

x=linspace(a,b,N); %待解x值

h=(b-a)/(N-1); %分成N-1个区间

y=zeros(1,N); %定义长度

y(1)=y0; %初始值

%具体的迭代过程

for i=1:N-1

k1=lkfun(x(i),y(i));

k2=lkfun(x(i)+h/2,y(i)+h*k1/2);

k3=lkfun(x(i)+h/2,y(i)+h*k2/2);

k4=lkfun(x(i)+h,y(i)+h*k3);

y(i+1)=y(i)+h/6*(k1+2*k2+2*k3+k4);

end

disp('常微分方程的数值解：')

disp([x;y])

lk(0,1,1,6); %调用已编写的龙格库塔函数LK

常微分方程的数值解：

0 0.2000 0.4000 0.6000 0.8000 1.0000

1.0000 0.8595 0.7676 0.7013 0.6505 0.6099

>> dsolve('Dy=-0.9*y/(1+2*t)','y(0)=1') %matlab自带求解函数dsolve，进行验证

ans =1/(2*t+1)^(9/20)

t=[0:0.2:1];

1./((2*t+1).^(9/20))

ans =

0 0.2000 0.4000 0.6000 0.8000 1.0000

1.0000 0.8595 0.7676 0.7013 0.6505 0.6100

%发现误差非常小

这是以前写的一个四阶的，自己修改一下应该可以用的