用四阶龙格-库塔方法求微分方程组

最近一段时间再忙期末考试，小学期课程设计的东西，没怎么更新博客....

更新一个用四阶龙格库塔方法求解脉冲微分方程，题目来源是一篇论文《Impulsive control of projective synchronization in chaotic systems 》。

概述：

做法其实也很简单，就是原来求一个因变量的龙格库塔公式变成求五个因变量的龙格库塔公式，当然，变化求解的过程是有不一样的地方的，需要特别处理一下。

% 1
t = [];
x = [];
y = [];
z = [];
sigma = 10;
u = 28;
p = 8 / 3;
X = [];
x(1) = 1;
y(1) = 2;
z(1) = 3;
X(1,1) = 2;
X(2,1) = 1;
h = 0.01;
a = 0;
b = 1;
t(1) = 0;
gain1 = -1.1;
gain2 = -1.1;
alpha = -3;
n = fix((b-a) / h);for i = 2 : n+1t(i) = t(i-1) +  h;k11 = sigma * ( y(i-1)  - x(i-1) );k21 = (u - z(i-1)) * x(i-1) - y(i-1);k31 = x(i-1) * y(i-1) - p * z(i-1);k41 = sigma * (  X(2, i-1) - X(1,i-1) );k51 = X(1,i-1) * (u - z(i-1)) - X(2,i-1);k12 = sigma * ( y(i-1) + h / 2 * k21 - x(i-1) - h / 2 * k11  );k22 = (u - z(i-1) - h / 2 * k31) * (x(i-1) + h / 2 * k11 ) - y(i-1) - h / 2 * k21;k32 = (x(i-1) + h / 2 * k11)* (y(i-1)+ h / 2 * k21 )- p * (z(i-1) + h /2 * k31  );k42 = sigma * (  X(2, i-1) + h/2*k51- X(1, i-1) - h/2*k41);k52 = (X(1,i-1) + h/2*k41  )* (u - z(i-1) - h/2*k31) - X(2,i-1) - h/2*k51;k13 = sigma * ( y(i-1) + h/2 * k22 -  x(i-1)  - h / 2 * k12 );k23 = (u - z(i-1) - h/2*k32 ) * (x(i-1)+h/2*k12) - y(i-1) - h / 2 * k22;k33 = (x(i-1) +h/2*k12) *( y(i-1)+ h/2 * k22 )- p *(z(i-1) + h / 2 * k32);k43 = sigma * (  X(2, i-1) + h/2*k52- X(1, i-1) - h/2*k42);k53 = (X(1,i-1) + h/2*k42  )* (u - z(i-1) - h/2*k32) - X(2,i-1) - h/2*k52;k14 = sigma * ( y(i-1)+h*k23 - x(i-1) - h * k13 );k24 = (u - z(i-1) - h*k33) * (x(i-1)+h*k13) - y(i-1) - h * k23;k34 = (x(i-1) +h*k13 ) * (y(i-1) + h*k23) - p * (z(i-1) + h * k33 );k44 = sigma * (  X(2, i-1) + h/2*k53- X(1, i-1) - h/2*k43);k54 = (X(1,i-1) + h/2*k43  )* (u - z(i-1) - h/2*k33) - X(2,i-1) - h/2*k53;x(i) = x(i-1) +  h / 6 * (k11 + 2 * k12 + 2 * k13 + k14); y(i) = y(i-1) + h / 6 * (k21 + 2 * k22 + 2 * k23 + k24);z(i) = z(i-1) + h / 6 * (k31 + 2 * k32 + 2 * k33 + k34);X(1, i) = X(1, i-1) + h /6 * (k41+ 2 * k42 + 2 * k42 + k44);X(2, i) = X(2, i-1) + h / 6 * (k51 + 2 * k52 + 2*k53 + k54);if mod(i, 15) == 0X(1, i) = X(1, i) + gain1 * (X(1,i) - alpha * x(i) );X(2, i) = X(2, i) + gain2 * (X(2,i) - alpha * y(i) );endend
%plot(t, x, 'b-', t, y, 'k-', t, z, 'c-', t, X(1,:), 'r-', t, X(2,:), 'm-')
% plot(t, x, t, y)
%  plot(t, X(1,:), t, X(2, :) )
e1 = X(1,:) - alpha * x;
e2 = X(2,:) - alpha * y;
plot(t, e1, t, e2)

matlab代码跑出来的图

用四阶龙格-库塔方法求微分方程组相关推荐

matlab：使用4阶龙格库塔方法求微分方程组的值
参考书目:<数值分析(matlab版)> 作者:周璐翻译 %调用龙格库塔方法求解微分方程组,P399,9.7 微分方程组 function [t,z] = my_rk4b(f, t0, ...
matlab：使用四阶龙格库塔方法求解微分方程组
%书籍:常用数值算法及其matlab实现 %第10章常微分方程初值问题的数值解法,例10.14使用 %四阶龙格库塔方法 function [t,z] = rk4symeq(fun, t0, tf, ...
计算方法实验（三）：四阶龙格-库塔方法
四阶Runge-Kutta数学原理给定常微分方程初值问题 {dydx=f(x,y),a≤x≤by(a)=αh=b−aN\left\{ \begin{matrix} \frac{\text{dy}}{ ...
matlab:使用4阶龙格库塔方法求解常微分方程组
%书籍:常用数值算法及其matlab实现 %第10章常微分方程初值问题的数值解法,例10.14使用 %四阶龙格库塔方法 function [t,z] = rk4symeq(fun, t0, tf, ...
[计算机数值分析]四阶龙格-库塔经典格式解常微分方程的初值问题
龙格-库塔方法的设计思想: 四阶龙格-库塔方法的经典格式: 程序设计框图: 例:设取步长h=0.2,从x=0到x=1用四阶经典格式解决以下常微分方程的初值问题. 运行示例: 源码: #include& ...
数值分析C++实现用四阶龙格-库塔(Runge-Kutta)方法求解常微分方程初值问题
问题:用四阶龙格-库塔(Runge-Kutta)方法求解常微分方程初值问题. 算法描述算法的程序框图: 具体算法: (1)读取a,b,n,f (2)计算步长h = (b-a)/n,x0=a,y0=f ...
四阶龙格库塔方程（Rungekutta）解二阶常微分方程组并计算船舶在迎浪下的纵摇埀荡耦合运动方程-附Matlab代码
今年年初的时候给师姐做了DDPG算法的船舶减横摇控制算法,师姐还有想法要让我把纵摇-埀荡两个自由度的减摇也做出来,这个任务归我了.实际上不管是多少个自由度的减摇,其实都需要解运动方程,当初做单自由度横 ...
matlab 数值计算课二阶微分方程-龙格库塔方法 ODE45
详见mathworks 龙格库塔方法写成矩阵(状态方程)的形式更简洁一点(其实这两种方法结果是一样的,如果C是[1,0,0]的话,就很明显了) 例如:求系统在0-5s内的单位阶跃相应,已知传递函数: ...
如何用龙格库塔方法求解范德波振子的运动
如何用龙格库塔方法求解范德波振子的运动坑,没有想好. 这是一个非常好的线索: Generalized Forced Van Der Pol Oscillator Phase Plot - fun ...

用四阶龙格-库塔方法求微分方程组

用四阶龙格-库塔方法求微分方程组相关推荐

最新文章

热门文章