sijie龙格库塔法matlab程序

1. 龙格库塔法matlab编程

用matlab的四、五阶龙格库塔方法（ode45）可以求得其数值解。设定其微分方程组的初始条件X(0)=0，S(0)=0。解决方法：

1、常微分方程组的自定义函数odefun（t，y），其主要内容

dy(1)=3*y(2)/(60+y(2))*y(1)-0.16*y(1);

dy(2)=4*250-4*y(2)-5*y(2)/(60+y(2))*y(1)-0.06*y(1);

2、t、y1、y2的求解

[t,y]=ode45(@odefun,[0 10],y0)

3、计算结果及t与S、X的关系图，S与X的关系图

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

%以下另存为文件myrk4.m
function[x,y]=myrk4(ufunc,y0,h,a,b)
%参数:函数名称，初始值向量，步长，时间起点，时间终点
n=floor((b-a)/h);%求步数
x(1)=a;%时间起点
y(:,1)=y0;%赋初值
%按龙格库塔方法进行求解
forii=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,[01],1);
subplot(212);plot(t,f);title('ode45求解结果')

3. 四阶常微分龙格—库塔法求解常微分方程的初值问题。 Matlab程序算例：用标准4级4阶R-K法求解

clear,clc	%清除内存中的变量
%数值解
y=inline('x*exp(x)+2*x-1');
y(1)

%四阶龙格库塔法
y0=[-132];
[x1,y1]=ode45(@fun,[0,1],y0);
y1(end,1)%y(1)的值

functiondy=fun(x,y)
dy=zeros(3,1);
dy(1)=y(2);
dy(2)=y(3);
dy(3)=y(3)+y(2)-y(1)+2*x-3;

运行结果

ans =

3.7183

ans =

3.7183