龙格库塔

最新推荐文章于 2024-10-13 12:14:17 发布

只会刷野的jug

最新推荐文章于 2024-10-13 12:14:17 发布

阅读量1.5k

点赞数 5

分类专栏： MATLAB 文章标签：龙格库塔最优控制1

MATLAB 专栏收录该内容

1 篇文章 1 订阅

订阅专栏

龙格库塔方法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);
  end
end

自己定义的微分方程组
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的曲线');```