龙格库塔方法MATLAB程序编写以及解决简单的微分方程组
%龙格库塔函数经典代码function [ x,y ] = Runge( dyfun,xspan,y0,h)x=xspan(1):h:xspan(2);y=zeros(length(y0),length(x));y(:,1) = y0(:); for n = 1:(length(x)-1) k1 = feval(dyfun,x(n),y(:,n)); k2 = feval(dyfun,x(n)+h/2,y(:,n)+h/2*k1); k3 = feval(dyfun,x(n)+h/2,y(:,n)+h/2*k2); k4 = feval(dyfun,x(n+1),y(:,n)+h*k3); y(:,n+1) = y(:,n) + h/6*(k1+2*k2+2*k3+k4); endend 自己定义的微分方程组function f = dyfun( t,y )a1 = 0.24;b1 = 0.21;b2 = 0.42;w = 0.15;f(1)=-a1*y(1)*(log(y(3)))/log(y(2)) - w*y(1)*log(y(2));f(2)= w*y(1)*log(y(2)) + b1*y(4)*log(y(2)) ;f(3)= b2*y(4);f(4)=a1*y(1)*(log(y(3)))/log(y(2)) - b1*y(4)*log(y(2)) - b2*y(4);f = f(:);end %编写注执行程序close all;clear all;clc;[x,y]=Runge(@dyfun,[0,10],[30,40,10,20],0.1);plot(x,y(1,:));title('控制前后S的曲线');figure;plot(x,y(2,:));title('控制前后A的曲线');figure;plot(x,y(3,:));title('控制前后I的曲线');figure;plot(x,y(4,:));title('控制前后C的曲线');```