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