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

龙格库塔法是一种求解高阶常微分方程的常用方法，在工程当中应用广泛，例如求解物体的运动方程等。
这里我们通过matlab程序编写龙格库塔算法求解二元常微分方程组，假设有常微分方程组：
{x¨−x˙+2y¨+y˙=−2sint−3costx¨+y¨=−sint−costx(0)=0y(0)=1x˙(0)=1y˙(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˙=−2sint−3costx¨+y¨=−sint−costx(0)=0y(0)=1x˙(0)=1y˙(0)=0
该方程有精确解：
{x=sinty=cost\left\{ \begin{array}{lr} x=\rm sin \it t \\ y=\rm cos \it t \\ \end{array} \right.{x=sinty=cost
令u1=xu_1=xu1=x,u2=x˙u_2=\dot{x}u2=x˙,w1=yw_1=yw1=y,w2=y˙w_2=\dot{y}w2=y˙,则原式可写为
{f1(t)=u2=x˙f2(t)=u2′=x¨=cost−x˙+y˙f3(t)=w2=y˙f4(t)=w2′=y¨=−sint−2cost+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¨=cost−x˙+y˙f3(t)=w2=y˙f4(t)=w2′=y¨=−sint−2cost+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)-1k11 = 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('计算值','精确值');

计算结果如图

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

matlab使用solve求解二元二次方程组
网上有些代码存在问题,这里只做了原理上的更正.也希望能各位码农能完全.准确地提供代码. 求解二元二次方程组代码: clc clear close all syms x y; f_1 = sym(x^2 ...
4 详解matlab实现龙格库塔算法求解复杂常微分方程组
目录 4.1 题目 4.2 问题背景 4.3 算法 4.4 matlab编程实现 4.5 结果展示 4.6 讨论 4.1 题目
Matlab求解常微分方程组
求解这个常微分方程组. 初始条件为其中ε取0.01,a是有上限的参数,求解方程的目的其实是找出a的临界值. syms y(t) for i = [0:0.5:1.5,1.7 ...
matlab龙格库塔法求通解,基于matlab及龙格库塔法求解布拉修斯方程.doc
基于matlab及龙格库塔法求解布拉修斯方程 Runge-Kutta法求解布拉修斯解摘要薄剪切层方程主要有三种解法,即相似解,非相似条件下对偏微分方程组的数值解和近似解.布拉修斯解是布拉修斯于19 ...
四阶龙格库塔法求解一次常微分方程组(python实现)
四阶龙格库塔法求解一次常微分方程组一.前言二.RK4求解方程组的要点 1. 将方程组转化为RK4求解要求的标准形式 2. 注意区分每个方程的独立性三.python实现RK4求解一次常微分方程组 ...
欧拉法、预估校正法(改进的欧拉法)与四阶龙格库塔法求解常微分方程的数值解C++程序
以y'=x+y,0<x<1,y(0)=1为例,取步长h=0.1,已知精确值为y=-x-1+2e^x,用来进行精度比较 #include<stdio.h> using names ...
python求解微分方程组_python – SymPy / SciPy：求解具有不同变量的常微分方程组...
我是SymPy和Python的新手,我目前正在使用Python 2.7和SymPy 0.7.5,其目标是: a)从文本文件中读取微分方程组 b)解决系统问题我已经阅读了this question和t ...
java常微分方程数值解,SymPy / SciPy：求解具有不同变量的常微分方程组
我是SymPy和Python的新手,我目前正在使用Python 2.7和SymPy 0.7.5,其目标是:a)从文本文件中读取微分方程组b)解决系统问题我已经读过this question和this ...
python求解微分方程组_python – SymPy / SciPy：求解具有不同变量的常微分方程组
我是SymPy和 Python的新手,我目前正在使用Python 2.7和SymPy 0.7.5,其目标是: a)从文本文件中读取微分方程组 b)解决系统问题我已经阅读了this question和 ...
一阶微分方程的物理意义_MIT—微分方程笔记24 一阶常微分方程组
18.03 微分方程 Differential Equations 第四单元一阶常微分方程组 Unit 4 First-Order Systems 第24讲一阶常微分方程组第25讲常系数齐次线 ...

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

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

最新文章

热门文章