MATLAB:微分方程(组)数值解

一、显式微分方程

 

clc,clear
tspan = [0:10]; 
y0 = 2;
[t1,y1] = ode23(@odefun_1,tspan,y0); %求数值解,精度相对低
[t2,y2] = ode113(@odefun_1,tspan,y0); %求数值解,精度相对高
yt = sqrt(tspan'+1)+1; %求精确解
subplot(1,2,1)
plot(t1,y1,'bo',t2,y2,'r*',tspan,yt,'g-')
legend('ode23求解','ode113求解','精确解','Location','best');
subplot(1,2,2)
y1t = abs(y1-yt); 
y2t = abs(y2-yt);
plot(t1,y1t,'bo',t2,y2t,'r*')
legend('ode23误差','ode113误差','Location','best');

function dydt = odefun_2(t,y)ft = 2*sin(t).*(t<4*pi) + 0.*(t>=4*pi); %定义分段函数f(t)gt = cos(t).*(t>=3.5*pi) + 0.*(t<3.5*pi); %定义分段函数g(t)dydt = [y(2) - ft;y(1)*gt-y(2)]; %微分方程组
end
ts = [0,20];
y0 = [1,2];
sol = ode23tb(@odefun_2,ts,y0)
subplot(1,2,1)
plot(sol.x,sol.y(1,:),'b-','LineWidth',1.5)
hold on
plot(sol.x,sol.y(2,:),'r--','LineWidth',1.5)
grid on
legend('{\ity}_1(t)','{\ity}_2(t)','Location','best')
legend('boxoff')
title('微分方程组解的曲线')
xlabel('t')
ylabel('y')
ysum = sum(sol.y);
subplot(1,2,2)
plot(sol.x,ysum,'b-','LineWidth',1.5)
hold on 
fplot(@(t)0*t,[0,20])
fh = @(t)sum(deval(sol,t));
[x01,y01] = fzero(fh,10)
[x02,y02] = fzero(fh,18)
plot(x01,y01,'r.','MarkerSize',15)
hold on
plot(x02,y02,'r.','MarkerSize',15)
grid on
legend('F(t) = {\ity}_1(t) + {\ity}_2(t)','y = 0的直线','Location','best');
legend('boxoff')

function dy = DyDtFun(t,y) dy = zeros(3,1);dy(1) = y(2)*y(3);dy(2) = -y(1)*y(3);dy(3) = -0.51*y(1)*y(2);
end
[T,Y]=ode45(@DyDtFun,[0 12],[0 1 1]);
plot(T,Y(:,1),'b:',T,Y(:,2),'g--',T,Y(:,3),'r.-','LineWidth',1.5)
grid on
legend('y_1','y_2','y_3')
legend('boxoff')
title('显示微分方程组初值问题求解')

 

function dydt = odefun_3(t,y) %定义函数文件dydt = [y(2);(1-y(1)^2)*y(2)-y(1)];
end
[t,y]=ode45(@odefun_3,[0 20],[2 0]);
plot(t,y(:,1),'b--',t,y(:,2),'r:','LineWidth',1.5);
grid on
title('常微分方程的解');
xlabel('t') 
ylabel('y')
legend('y(t)','dy/dt','Location','best')
legend('boxoff')

 

function dxdt = apolloeq(t,x)mu = 1/82.45;mu1 = 1 - mu;r1 = sqrt((x(1)+mu).^2 + x(3).^2);r2 = sqrt((x(1)-mu1).^2 + x(3).^2);dxdt = [x(2);2*x(4) + x(1) - mu1*(x(1) + mu)/r1.^3 - mu*(x(1) - mu1)/r2.^3;x(4);-2*x(2) + x(3) - mu1*x(3)/r1.^3 - mu*x(3)/r2.^3];
end
x0 = [1.2;0;0;-1.04935751];
% [t,y] = ode45(@apolloeq,[0,20],x0);
[t,y] = ode45(@apolloeq,[0,20],x0,odeset('RelTol',1e-6));
plot(y(:,1), y(:,3),'LineWidth',2)
grid on
title('Apollo卫星的运动轨迹求解图像')

看课后案例分析。 

二、完全隐式微分方程

function dydt = odeifun_1(t,y,yp)dydt = t.*y.^2.*yp.^3 - y.^3.*yp.^2 + t.*(t.^2 + 1).*yp - t.^2.*y;
end
clc,clear
t0 = 1; 
y0 = sqrt(3/2); 
ypo = 0;
[y0,yp0] = decic(@odeifun_1,t0,y0,1,ypo,0)
%求得y0 = 1.2247 yp0 = 0.8165
[t,y] = ode15i(@odeifun_1,[1,10],y0,yp0,odeset('RelTol',1e-5));
plot(t,y)
hold on
plot(t,sqrt(t.^2+0.5),'r*') %绘制精确解散点
legend('y=ode15i','y=sqrt(t.^2+0.5)')

function dydt = odeifun_2(t,y,dy)
dydt = [dy(1)-y(2);dy(2)*sin(y(4))+dy(4)^2 + 2*y(1)*y(3)-y(1)*dy(2)*y(4);dy(3)-y(4);y(1)*dy(2)*dy(4)+cos(dy(4))-3*y(2)*y(3)];
end
t0 = 0; 
y0 = [1;0;0;1]; 
fix_y0 = ones(4,1); 
dy0 = [0 3 1 0]';
fix_dy0 = zeros(4,1); 
[y02,dy02] = decic(@odeifun_2,t0,y0,fix_y0,dy0,fix_dy0); [t y] = ode15i(@odeifun_2,[0 5],y02,dy02); 
plot(t,y(:,1),'k-',t,y(:,2),'r--', t,y(:,3),'g-.',t,y(:,4),'b:','linewidth',2)
grid on
legend('{\ity}_1(t)','{\ity}_2(t)','{\ity}_3(t)','{\ity}_4(t)','Location','best')
title('隐式微分方程组的解曲线')
xlabel('t')
ylabel('y')
legend('boxoff')

三、代数微分方程组(DAEs) 

ode15s用于求显式的微分方程组

y0 = [0.8;0.1;0.1]; %初值
M = [1 0 0;0 1 0;0 0 0]; %质量矩阵
options = odeset('Mass',M,'RelTol',1e-4,'AbsTol',[1e-6,1e-6,1e-6]);
tspan = [0 20];
[T,Y] = ode15s(@DAEsfun_1,tspan,y0,options);
plot(T,Y(:,1),'b-','linewidth',2)
hold on
plot(T,Y(:,2),'r-.','linewidth',2)
plot(T,Y(:,3),'g:','linewidth',2)
legend('{\ity}_1(t)','{\ity}_2(t)','{\ity}_3(t)','Location','best');
title('代数微分方程求解结果图')

ode15i用于求隐式的微分方程组 

function dydt = iDAEsfun_2(t,y,yp)% 注意与显式微分方程组定义的不同dydt = [yp(1)+0.3*y(1)+2*y(2).*sin(yp(3))+y(2).*y(4);yp(2)+y(2)+0.5*cos(yp(1)+y(3))+0.2*sin(0.6*t);yp(3)+0.2*y(1).*y(2)-exp(-yp(1));y(1)+y(2)-y(3)-y(4)-1];
end
y0 = [1;0.5;0.3;0.2]; % 初值
tspan = [0 5]; % 求值区间
fix_y0 = [0;1;1;1];
fix_yp0 = [0;0;0;0]; % 该组初值都可以改变,故全部为0
% fsolve
fh = @(y)[y(1)+0.4+sin(y(3)),y(2)+0.5+0.5*cos(y(1)+0.3),y(3)+0.1-exp(y(1))];
fsolve(fh,[1,1,1]) % -0.7596 -0.9481  0.3678
yp0 = [-0.8;-1;1;1];
% 首先求yp0,使得f(t,y,yp) = 0
[y02,yp02] = decic(@iDAEsfun_2,0,y0,fix_y0,yp0,fix_yp0)
[t,y] = ode15i(@iDAEsfun_2,tspan,y02,yp02);
plot(t,y(:,1),'k-',t,y(:,2),'r-.',t,y(:,3),'b:',t,y(:,4),'g--','linewidth',1.5)
grid on
legend('{\ity}_1(t)','{\ity}_2(t)','{\ity}_3(t)','{\ity}_4(t)','Location','best')
legend('boxoff')

四、延时/时滞微分方程 

 

显式微分方程 

function dydt = DDEsfun_1(t,y,Z)dydt = [0.5*Z(3,2) + 0.5*y(2).*cos(t);0.3*Z(1,1) + 0.7*y(3).*sin(t);y(2) + cos(2*t)];
end
lags = [1,3]; % 延迟常数向量
history = [0,0,1]; % 小于初值时的历史函数
tspan = [0,8];
sol = dde23(@DDEsfun_1,lags,history,tspan) % 调用dde23求解
ti = 0:0.01:8;
% yi = deval(sol,ti);
% plot(ti,yi(1,:),ti,yi(2,:),ti,yi(3,:))
plot(sol.x,sol.y(1,:),'r-','linewidth',2);
hold on
grid on
plot(sol.x,sol.y(2,:),'b-.','linewidth',2)
plot(sol.x,sol.y(3,:),'g-*','linewidth',1)
legend('{\ity}_1(t)','{\ity}_2(t)','{\ity}_3(t)','Location','best')
title('方程各解的曲线图')
legend('boxoff')

function dydt = DDEsfun_2(t,x,Z)dydt = [-2*x(2) - 3*Z(1,1);-0.05*x(1).*x(3) - 2*Z(2,2) + 2;0.3*x(1).*x(2).*x(3) + cos(x(1).*x(2)) + 2*sin(0.1*t.^2)];
end
lags = @(t,x)[t - 0.2*abs(sin(t)); t - 0.8];
history = zeros(3,1);
tspan = [0,10];
sol = ddesd(@DDEsfun_2,lags,history,tspan)
plot(sol.x,sol.y(1,:),'r--',sol.x,sol.y(2,:),'b-.',sol.x,sol.y(3,:),'g:','linewidth',2);
legend('{\ity}_1(t)','{\ity}_2(t)','{\ity}_3(t)','Location','best')
legend('boxoff')
grid on
title('变时间延迟微分方程各解的曲线图')

lags = @(t,x)[t - 0.2*abs(sin(t)); t - 0.8];
% history = zeros(3,1);
history = @(t,x)[sin(t+1); cos(t); exp(3*t)];
tspan = [0,10];
sol = ddesd(@DDEsfun_2,lags,history,tspan)
% ti = 0:0.01:10;
% yi = deval(sol,ti);
% plot(ti,yi(1,:),'r--',ti,yi(2,:),'b-.',ti,yi(3,:),'g:','linewidth',2)
plot(sol.x,sol.y(1,:),'r--',sol.x,sol.y(2,:),'b-.',sol.x,sol.y(3,:),'g:','linewidth',2)
legend('{\ity}_1(t)','{\ity}_2(t)','{\ity}_3(t)','Location','best')
legend('boxoff')
grid on
title('变时间延迟微分方程各解的曲线图')

function dxdt = DDEsfun_3(t,x,Z,ZP)A1 = [-13 3 -3; 106 -116 62; 207 -207 113];A2 = diag([0.02 0.03 0.04]);B = [0;1;2];U = 1;dxdt = A1*Z + A2*ZP + B*U;
end
history = zeros(3,1);
sol = ddensd(@DDEsfun_3,0.15,0.5,history,[0,15])
plot(sol.x,sol.y(1,:),'r--',sol.x,sol.y(2,:),'b-.',sol.x,sol.y(3,:),'g:','linewidth',2)
legend('{\ity}_1(t)','{\ity}_2(t)','{\ity}_3(t)','Location','best')
legend('boxoff')
grid on
title('中立型延迟微分方程组解曲线')

五、微分方程边值问题

 solinit: solve initial.

 

function dydt = bvpfun_1(t,y)dydt = [y(2);2*y(2).*cos(t) - y(1).*sin(4*t) - cos(3*t)];
end
function res = bcfun_1(ya,yb)res = [ya(1) - 1;yb(1) - 2];
end
function yint = initfun_1(t)%  y(0) = 1, y(4) = 2 --> (0,1),(4,2)yint = [1/4*t + 1;1/4]
end
function dfdy = jacb(t,y)dfdy = [0, 1;-sin(4*t), 2*cos(t)];
end
solinit = bvpinit(linspace(0,4,10),@initfun_1)
opts = bvpset('FJacobian',@jacb,'RelTol',1e-4,'AbsTol',1e-8,'Stats','on');
sol = bvp4c(@bvpfun_1,@bcfun_1,solinit,opts)
sol5 = bvp5c(@bvpfun_1,@bcfun_1,solinit,opts)
% bvp5c精度比bvp4c高
ti = 0:0.01:4;
yi = deval(sol,ti);
plot(ti,yi(1,:),ti,yi(2,:))
grid on
legend('{\ity}_1(t)','{\ity}_2(t)','Location','best')
legend('boxoff')
title('常微分方程边界问题解的曲线图')

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.mzph.cn/news/773675.shtml

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

C语言:动态内存管理(malloc,calloc,realloc,free)

目录 前言 malloc函数 free函数 calloc函数 realloc函数 前言 在这一章节将讲解动态内存分配&#xff0c;它可以在程序的堆区创建一块内存&#xff0c;在这块内存中存什么值就是由自己决定的了 开辟的空间有两个特点&#xff1a; 1. 空间开辟的大小是固定的 2. 数组在…

视频 | 轨迹模型及其它潜变量模型理论与实践

2024年3月20日&#xff0c;郑老师开办了“春分”临床统计学沙龙&#xff0c;本期沙龙和诸位分享了轨迹增长模型&#xff0c;以及其它潜变量模型包括潜类别模型和潜剖面模型的详细介绍和文献解读&#xff0c;R语言代码实践也有分享&#xff01; 非常感谢各位朋友在忙碌的工作生活…

线性数据结构----(数组,链表,栈,队列,哈希表)

线性数据结构 数组链表栈使用场景 队列应用场景 哈希表特点哈希函数&#xff0c;哈希值&#xff0c;哈希冲突键值对 Entry 开放寻址法和拉链法 参考文档 数组 数组(Array) 是一种很常见的数据结构。由相同类型的元素组成&#xff0c;并且是使用一块连续的内存来存储的。 在数组…

安装redis时候修改过的配置文件

只要是石头&#xff0c;到哪里都不会发光的 bind 绑定主机某个网卡对应的IP地址&#xff0c;如果某个主机有两个网卡A和B&#xff0c;那么绑定了A&#xff0c;通过B连接就会无法访问protected-mode 保护模式 Yes为只能本地访问port 启动的端口号pidfile pid存放的位置&#xff…

python django实战开发序列化器的一个应用心得分享

需求: 查询的时候返回不包括SharePasswd 字段, 但是新增操作需要用到该字段 再不写多个model模型和序列化器的前提下实现 如果您在查询&#xff08;GET 请求&#xff09;时不希望返回 SharePasswd 字段&#xff0c;但在新增&#xff08;POST 请求&#xff09;时需要用到该字段…

Python3爬取2023省市区

爬取地址https://www.stats.gov.cn/sj/tjbz/tjyqhdmhcxhfdm/2023/ import re import requests import pandas as pd import warnings warnings.filterwarnings("ignore") import time from lxml import etree import pymysql t ,urls ,names [],[],[] INDEX_URL &…

Java两地经纬度通过高德api获取两地距离(公里)

代码如下&#xff1a; String startLongitude entity.getLONGITUDE(); // 起点&#xff08;当前位置&#xff09;经度String startLatitude entity.getLATITUDE(); // 起点纬度String endLongitude entity.getLO(); // 终点经度String endLatitude entity.getLA(); …

Spring框架介绍及详细使用

前言 本篇文章将会对spring框架做出一个比较详细的讲解&#xff0c;并且每个知识点基本都会有例子演示&#xff0c;详细记录下了我在学习Spring时所了解到全部知识点。 在了解是什么spring之前&#xff0c;我们要先知道spring框架在开发时&#xff0c;服务器端采用三层架构的方…

电子元器件批发采购的数字化转型与技术应用

电子元器件批发采购的数字化转型和技术应用是提高效率、优化供应链管理以及增强竞争力的关键。以下是一些数字化转型和技术应用的方面&#xff1a; 电子商务平台&#xff1a;建立并优化电子商务平台&#xff0c;提供在线采购功能&#xff0c;使客户能够方便地浏览产品、下订单、…

ABNDP: Co-optimizing Data Access and Load Balance in Near-Data Processing——论文泛读

ASPLOS 2023 Paper 论文阅读笔记整理 问题 近数据处理&#xff08;NDP&#xff09;是一种很有前途的体系结构范式&#xff0c;可以解决数据密集型应用程序的内存墙挑战。基于3D堆叠存储器的典型NDP系统包含大量并行处理单元&#xff0c;每个并行处理单元都可以访问其本地存储…

HTML基础:8个常见表单元素的详解

你好&#xff0c;我是云桃桃。 一个希望帮助更多朋友快速入门 WEB 前端程序媛。 后台回复“前端工具”可免费获取开发工具&#xff0c;持续更新。 今天来说说 HTML 表单。它是用于收集用户输入信息的元素集合。例如文本框、单选按钮、复选框、下拉列表等。 用户经常填写的表…

​深入理解JMeter性能测试日志:分析并发用户行为与吞吐量指标

在进行软件性能测试时&#xff0c;使用Apache JMeter等工具能够生成详尽的测试报告和日志。这些数据对于评估应用程序的性能至关重要。本文将针对JMeter生成的压测日志数据进行详细解析&#xff0c;并解释特定字段的意义。 日志数据摘要&#xff1a; 以下是整个JMeter压测日志…

2024智能EDM邮件营销系统使用攻略

在数字化营销领域&#xff0c;智能EDM&#xff08;Electronic Direct Mail&#xff09;邮件营销作为一种高效、精准的推广方式&#xff0c;正日益受到企业的高度重视。而要实现这一策略的成功落地&#xff0c;一个高可靠性和高稳定性的专业邮件发送平台则是不可或缺的关键环节。…

前后端实时数据通信

实现前后端实时数据转换通常涉及到以下几个步骤&#xff1a; 后端提供数据转换接口。 前端实时数据获取。 前端实时数据转换。 前端实时展示转换后数据。 以下是一个简单的例子&#xff0c;假设后端提供了一个接口来转换某种数据格式&#xff0c;前端使用JavaScript和WebS…

大数据分析案例-基于决策树算法构建大学毕业生薪资预测模型

&#x1f935;‍♂️ 个人主页&#xff1a;艾派森的个人主页 ✍&#x1f3fb;作者简介&#xff1a;Python学习者 &#x1f40b; 希望大家多多支持&#xff0c;我们一起进步&#xff01;&#x1f604; 如果文章对你有帮助的话&#xff0c; 欢迎评论 &#x1f4ac;点赞&#x1f4…

前端发版上线出现白屏问题

目录 路由配置问题资源缓存问题首屏加载过慢 &#xff1a;喂&#xff0c;你的页面白啦&#xff01; 出现上线白屏的问题有很多&#xff0c;如&#xff1a;配置错误、缓存问题、浏览器兼容问题&#xff0c;根据不同情况去解决。 路由配置问题 问题描述&#xff1a; 在vue开发…

C语言中位运算介绍

在C语言中&#xff0c;位运算是一种对二进制位进行操作的运算方式&#xff0c;它可以对数据的二进制表示进行位级别的操作&#xff0c;包括按位与、按位或、按位异或、按位取反等。位运算常用于处理底层数据结构、优化代码性能以及实现各种算法。本文将深入介绍C语言中的位运算…

两区域二次调频风火机组,麻雀启发式算法改进simulink与matlab联合

区域1结果 区域2结果 红色曲线为优化后结果〔风火机组二次调频〕

软件杯 深度学习+opencv+python实现车道线检测 - 自动驾驶

文章目录 0 前言1 课题背景2 实现效果3 卷积神经网络3.1卷积层3.2 池化层3.3 激活函数&#xff1a;3.4 全连接层3.5 使用tensorflow中keras模块实现卷积神经网络 4 YOLOV56 数据集处理7 模型训练8 最后 0 前言 &#x1f525; 优质竞赛项目系列&#xff0c;今天要分享的是 &am…

鸿蒙操作系统-初识

HarmonyOS-初识 简述安装配置hello world1.创建项目2.目录解释3.构建页面4.真机运行 应用程序包共享包HARHSP 快速修复包 官方文档请参考&#xff1a;HarmonyOS 简述 1.定义&#xff1a;HarmonyOS是分布式操作系统&#xff0c;它旨在为不同类型的智能设备提供统一的操作系统&a…