Matlab 四阶龙格库塔法求解二元常微分方程组

常用小程序 同时被 2 个专栏收录
3 篇文章 0 订阅
6 篇文章 0 订阅

龙格库塔法是一种求解高阶常微分方程的常用方法,在工程当中应用广泛,例如求解物体的运动方程等。
这里我们通过matlab程序编写龙格库塔算法求解二元常微分方程组,假设有常微分方程组:
{ x ¨ − x ˙ + 2 y ¨ + y ˙ = − 2 s i n t − 3 c o s t x ¨ + y ¨ = − s i n t − c o s t x ( 0 ) = 0 y ( 0 ) = 1 x ˙ ( 0 ) = 1 y ˙ ( 0 ) = 0 \left\{ \begin{array}{lr} \ddot{x}-\dot{x}+2\ddot{y}+\dot{y}=-2\rm sin \it t - \rm 3\rm cos \it t \\ \ddot{x}+\ddot{y}=-\rm sin \it t - \rm \rm cos \it t \\ x(0)=0 \\ y(0)=1 \\ \dot{x}(0)=1 \\ \dot{y}(0)=0 \\ \end{array} \right. x¨x˙+2y¨+y˙=2sint3costx¨+y¨=sintcostx(0)=0y(0)=1x˙(0)=1y˙(0)=0
该方程有精确解:
{ x = s i n t y = c o s t \left\{ \begin{array}{lr} x=\rm sin \it t \\ y=\rm cos \it t \\ \end{array} \right. {x=sinty=cost
u 1 = x u_1=x u1=x, u 2 = x ˙ u_2=\dot{x} u2=x˙, w 1 = y w_1=y w1=y, w 2 = y ˙ w_2=\dot{y} w2=y˙,则原式可写为
{ f 1 ( t ) = u 2 = x ˙ f 2 ( t ) = u 2 ′ = x ¨ = c o s t − x ˙ + y ˙ f 3 ( t ) = w 2 = y ˙ f 4 ( t ) = w 2 ′ = y ¨ = − s i n t − 2 c o s t + x ˙ − y ˙ \left\{ \begin{array}{lr} f1(t)=u_2=\dot{x} \\ f2(t)=u_2'=\ddot{x}=\rm cos \it t -\dot{x}+\dot{y} \\ f3(t)=w_2=\dot{y} \\ f4(t)=w_2'=\ddot{y}=-\rm sin \it t - \rm 2\rm cos \it t +\dot{x}-\dot{y} \\ \end{array} \right. f1(t)=u2=x˙f2(t)=u2=x¨=costx˙+y˙f3(t)=w2=y˙f4(t)=w2=y¨=sint2cost+x˙y˙
可以编写四个子函数如下:

function output = f1(x,u1,u2,w1,w2)
output = u2;
function output = f2(x,u1,u2,w1,w2)
output = cos(x)-u2+w2;
function output = f3(x,u1,u2,w1,w2)
output = w2;
function output = f4(x,u1,u2,w1,w2)
output = -sin(x)-2*cos(x)-w2+u2;

四阶龙格库塔函数程序

function [u1,u2,w1,w2] = RK4_2variable(u1,u2,w1,w2,h,a,b)

x = a:h:b;

for i = 1:length(x)-1

k11 = f1(x(i) , u1(i) , u2(i) , w1(i) , w2(i));
k21 = f2(x(i) , u1(i) , u2(i) , w1(i) , w2(i));
L11 = f3(x(i) , u1(i) , u2(i) , w1(i) , w2(i));
L21 = f4(x(i) , u1(i) , u2(i) , w1(i) , w2(i));

k12 = f1(x(i)+h/2 , u1(i)+h*k11/2 , u2(i)+h*k21/2 , w1(i)+h*L11/2, w2(i)+h*L21/2);
k22 = f2(x(i)+h/2 , u1(i)+h*k11/2 , u2(i)+h*k21/2 , w1(i)+h*L11/2, w2(i)+h*L21/2);
L12 = f3(x(i)+h/2 , u1(i)+h*k11/2 , u2(i)+h*k21/2 , w1(i)+h*L11/2, w2(i)+h*L21/2);
L22 = f4(x(i)+h/2 , u1(i)+h*k11/2 , u2(i)+h*k21/2 , w1(i)+h*L11/2, w2(i)+h*L21/2);

k13 = f1(x(i)+h/2 , u1(i)+h*k12/2 , u2(i)+h*k22/2 , w1(i)+h*L12/2 , w2(i)+h*L22/2);
k23 = f2(x(i)+h/2 , u1(i)+h*k12/2 , u2(i)+h*k22/2 , w1(i)+h*L12/2 , w2(i)+h*L22/2);
L13 = f3(x(i)+h/2 , u1(i)+h*k12/2 , u2(i)+h*k22/2 , w1(i)+h*L12/2 , w2(i)+h*L22/2);
L23 = f4(x(i)+h/2 , u1(i)+h*k12/2 , u2(i)+h*k22/2 , w1(i)+h*L12/2 , w2(i)+h*L22/2);

k14 = f1(x(i)+h , u1(i)+h*k13 , u2(i)+h*k23 , w1(i)+h*L13 , w2(i)+h*L23);
k24 = f2(x(i)+h , u1(i)+h*k13 , u2(i)+h*k23 , w1(i)+h*L13 , w2(i)+h*L23);
L14 = f3(x(i)+h , u1(i)+h*k13 , u2(i)+h*k23 , w1(i)+h*L13 , w2(i)+h*L23);
L24 = f4(x(i)+h , u1(i)+h*k13 , u2(i)+h*k23 , w1(i)+h*L13 , w2(i)+h*L23);

u1(i+1) = u1(i) + h/6 * (k11 + 2*k12 + 2*k13 + k14);
u2(i+1) = u2(i) + h/6 * (k21 + 2*k22 + 2*k23 + k24);
w1(i+1) = w1(i) + h/6 * (L11 + 2*L12 + 2*L13 + L14);
w2(i+1) = w2(i) + h/6 * (L21 + 2*L22 + 2*L23 + L24);

end
end

计算主程序

% main func
clear;clc;
u1(1) = 0;
u2(1) = 1;
w1(1) = 1;
w2(1) = 0;
h=0.01;
a = 0;b=20;
[u1,u2,w1,w2] = RK4_2variable(u1,u2,w1,w2,h,a,b);
figure 
plot(a:h:b,u1,'r-');
hold on
plot(a:h:b,sin(a:h:b),'b-.');
xlabel('time');
ylabel('x');
legend('计算值','精确值');

计算结果如图
在这里插入图片描述

  • 8
    点赞
  • 2
    评论
  • 31
    收藏
  • 一键三连
    一键三连
  • 扫一扫,分享海报

相关推荐
©️2020 CSDN 皮肤主题: 大白 设计师:CSDN官方博客 返回首页
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、C币套餐、付费专栏及课程。

余额充值