matlab:使用四阶龙格库塔方法求解微分方程组
%书籍:常用数值算法及其matlab实现
%第10章 常微分方程初值问题的数值解法,例10.14使用
%四阶龙格库塔方法
function [t,z] = rk4symeq(fun, t0, tf, Za, h)
%fun:微分方程的右表达式
%t0, tn为区间
%Za为初值,是列向量
M = floor(tf-t0)/h ; %离散点的个数M+1
if t0 >= tfprintf('左端点必须小于右端点');return;
end
N = length(Za); %获得变量个数,N
z = zeros(M+1, N);
t = zeros(M +1, 1);
t =[t0 : h :tf]';
z(1,:) = Za'; %假设Za为列向量,与微分方程中的变量方向统一,变成行向量for i = 1:MK1 = feval(fun, t(i) , z(i,:)); %K是行向量K2 = feval(fun, t(i)+1/2*h ,z(i,:)+1/2* h*K1);K3 = feval(fun, t(i)+1/2*h ,z(i,:)+1/2* h*K2);K4 = feval(fun, t(i)+ h ,z(i,:)+ h*K3); z(i+1,:) = z(i,:) +h/6 *(K1 + 2*K2 + 2*K3 + K4);
end
方程组1-输出为行向量
%书籍:常用数值算法及其matlab实现
%第10章 常微分方程初值问题的数值解法
%四阶龙格库塔,例10.15
function s = exa10_15(t,z)
%z是个向量,1*3
%输出s也是向量,1*3
s = zeros(1,3);
dy1 = -10*z(1) +10*z(2);
dy2 = 28*z(1) - z(2) - z(1)*z(3);
dy3 = -8/3*z(3) + z(1)*z(2);
s = [dy1 dy2 dy3]; %输出s为行向量微分方程组——输出为列向量
%书籍:常用数值算法及其matlab实现
%第10章 常微分方程初值问题的数值解法
%四阶龙格库塔,例10.15
function s = exa10_15b(t,z)
%z是个向量,1*3
%输出s也是向量,1*3
s = zeros(1,3);
dy1 = -10*z(1) +10*z(2);
dy2 = 28*z(1) - z(2) - z(1)*z(3);
dy3 = -8/3*z(3) + z(1)*z(2);
s = [dy1; dy2; dy3]; %输出s为列向量,用于ode45调用
主函数
clear all;clc;close all;
%书籍:常用数值算法及其matlab实现
%第10章 常微分方程组初值问题的数值解法
%四阶龙格库塔,例10.15
%函数原型 function [t,z] = rk4symeq(fun, t0, tn, Za, h)
format long
t0 = 0; tf = 5;
Za = [-8 ; 8; 27]; %x初值
h1 = 0.01; [t1,z1] = rk4symeq(@exa10_15, t0, tf , Za, h1);
figure(1)
plot(t1,z1(:,1),'b',t1,z1(:,2), 'r',t1,z1(:,3), 'g--')
legend('y1','y2','y3')
figure(2)
plot3(z1(:,1),z1(:,2),z1(:,3));
xlabel('x');ylabel('y');zlabel('z');[t2,z2] = ode45(@exa10_15b, [0,5] ,[-8 , 8 ,27]); %用系统函数
figure(3)
plot(t2,z2(:,1),'b',t2,z2(:,2), 'r',t2,z2(:,3), 'g--')
legend('y1','y2','y3')
运行结果:
matlab:使用四阶龙格库塔方法求解微分方程组相关推荐
- 用四阶龙格-库塔方法求微分方程组
最近一段时间再忙期末考试,小学期课程设计的东西,没怎么更新博客.... 更新一个用四阶龙格库塔方法求解脉冲微分方程,题目来源是一篇论文<Impulsive control of projecti ...
- matlab:使用4阶龙格库塔方法求微分方程组的值
参考书目:<数值分析(matlab版)> 作者:周璐 翻译 %调用龙格库塔方法求解微分方程组,P399,9.7 微分方程组 function [t,z] = my_rk4b(f, t0, ...
- matlab:使用4阶龙格库塔方法求解常微分方程组
%书籍:常用数值算法及其matlab实现 %第10章 常微分方程初值问题的数值解法,例10.14使用 %四阶龙格库塔方法 function [t,z] = rk4symeq(fun, t0, tf, ...
- 如何用龙格库塔 方法求解 范德波振子的运动
如何用龙格库塔 方法求解 范德波振子的运动 坑,没有想好. 这是一个非常好的线索: Generalized Forced Van Der Pol Oscillator Phase Plot - fun ...
- 计算方法实验(三):四阶龙格-库塔方法
四阶Runge-Kutta数学原理 给定常微分方程初值问题 {dydx=f(x,y),a≤x≤by(a)=αh=b−aN\left\{ \begin{matrix} \frac{\text{dy}}{ ...
- 欧拉方法求解微分方程C语言,用欧拉方法求解微分方程组
将来,如果您在问题中包含完整的警告消息,这将很有帮助-它将包含发生问题的行:tmp/untitled.py:15: RuntimeWarning: overflow encountered in do ...
- matlab 龙格-库塔 法求解常微分方程
最近学习分室模型,里面碰到了用matlab 龙格-库塔 法求解常微分方程 研究了一阵子终于明白到底怎么实现了: 1. matlab 新建.m文件,编写龙格-库塔法求解函数 function [x,y] ...
- 【Hydro】龙格-库塔方法的公式推导
一.龙格-库塔方法的公式推导 摘自<数值计算方法(MATLAB版)> 二.龙格-库塔方法在水库调洪演算的应用 一般库水面坡降很小,忽略动库容影响,近似看成静水面,水库蓄水量V只随坝前水位Z ...
- matlab:使用龙格库塔法求解微分方程组
%书籍:常用数值算法及其matlab实现 %第10章 常微分方程初值问题的数值解法,例10.14使用 %四阶龙格库塔方法 function [t,z] = rk4symeq(fun, t0, tf, ...
最新文章
- Percona XtraBackup备份到恢复记录
- 爬虫python下载-python 下载 爬虫
- Vue---淘宝镜像cnpm---安装vue-cli
- xyntservice
- linux磁盘空间回收,linux磁盘空间满了后怎么去判定哪个地方占了多大的空间并回收...
- jQuery Mobile中按钮<a>或<button>添加class样式ui-btn-*
- CentOS7 配置免密登陆
- VS2013 工具包修复记录
- ofd 文件发票解析
- 华为设备DHCP/DHCP中继及DHCPv6/DHCPv6中继详解
- 后台 事务SqlTransaction
- 使用bootloader进行远程固件升级(32MCU)
- PC-DMIS 圆(环)型平面自动测量子程序
- 最简单的正版Win10重装系统教程
- xp系统 共享文件夹服务器,WinXP如何共享文件夹?共享文件夹的方法
- 03 循环、字符串、列表、元祖、字典
- 对UART、RS232、485通信的理解
- ZOJ1084(四色定理)
- Pspice17.2绘制原理图
- hdu 5107 线段树+离散化+归并排序+极角排序