Question

急求MATLAB编程源代码用四阶龙格库塔法解如下微分方程 y'=y-2x/y(0<x<1),y(0)=1,步长为h=0.2

不允许调用ode45函数，要求自己编写龙格库塔过程

Answer 1

% 以下另存为文件 myrk4.m
function [x,y]=myrk4(ufunc,y0,h,a,b)
%参数: 函数名称，初始值向量，步长，时间起点，时间终点
n=floor((b-a)/h);%求步数 
x(1)=a;%时间起点  
y(:,1)=y0;%赋初值
%按龙格库塔方法进行求解 
for ii=1:n   
	x(ii+1)=x(ii)+h;   
	k1=ufunc(x(ii),y(:,ii));   
	k2=ufunc(x(ii)+h/2,y(:,ii)+h*k1/2);   
	k3=ufunc(x(ii)+h/2,y(:,ii)+h*k2/2); 
	k4=ufunc(x(ii)+h,y(:,ii)+h*k3);
	y(:,ii+1)=y(:,ii)+h*(k1+2*k2+2*k3+k4)/6; 
end

以下是主程序

% y'=y-2x/y  (0<x<1),y(0)=1,步长为h=0.2
fun = inline('y-2*x/y');
[t1,f1]=myrk4(fun,1,0.2,0,1);%测试时改变test_fun的函数维数，别忘记改变初始值的维数 
subplot(211); plot(t1,f1)  %自编函数 
title('自编函数求解结果') 
%用系统自带函数ode45进行比较
[t,f] = ode45(fun,[0 1],1); 
subplot(212); plot(t,f);title('ode45求解结果')