【数学建模】微分方程的数值求解

微分方程的数值求解

  • 一阶差分求解微分方程原理:
  • 四阶龙格-库塔方法
    • 应用:小船渡河问题:
  • 进阶求二阶微分方程

一阶差分求解微分方程原理:

d y d x = f ( x n , y n ) \dfrac{dy}{dx}=f(x_n,y_n) dxdy=f(xn,yn)

y n + 1 − y n x n + 1 − x n = f ( x n , y n ) \dfrac{y_{n+1}-y_n}{x_{n+1}-x_n}=f(x_n,y_n) xn+1xnyn+1yn=f(xn,yn)

h = x n + 1 − x n h = x_{n+1} - x_n h=xn+1xn

y n + 1 = y n + h f ( x n , y n ) y_{n+1}=y_n+hf(x_n,y_n) yn+1=yn+hf(xn,yn)

例子:
{ d y d x = x y y 0 = 1 \begin{cases} \dfrac{dy}{dx}=xy\\ y_0 = 1 \end{cases} dxdy=xyy0=1

通过高等数学可以求得解析解为: y = e 1 2 x 2 y=e^{\frac{1}{2}x^2} y=e21x2

如果通过一阶差分求解
{ y n + 1 = y n + h x n y n x n + 1 = x n + h y 0 = 1 x 0 = 0 \begin{cases} y_{n+1}=y_n+hx_ny_n\\ x_{n+1}=x_n+h\\ y_0 = 1\\ x_0 = 0 \end{cases} yn+1=yn+hxnynxn+1=xn+hy0=1x0=0

matlab求解:

y = zeros(1,1000);
x = zeros(1,1000);
y(1)=1;
x(1)=0;
h=0.1;for n=1:1000y(n+1)=y(n)+h*x(n)*y(n);x(n+1)=x(n)+h;
end
plot(x,y,'-');
xlim([0,6]);
ylim([0,100]);

在这里插入图片描述
如果再加上解析解

y = zeros(1,1000);
x = zeros(1,1000);
y(1)=1;
x(1)=0;
h=0.1;for n=1:1000y(n+1)=y(n)+h*x(n)*y(n);x(n+1)=x(n)+h;
end
plot(x,y,'-');
hold on;
ezplot('exp(x^2/2)');
xlim([0,6]);
ylim([0,100]);

在这里插入图片描述
h=0.01如下图
在这里插入图片描述

四阶龙格-库塔方法

原理/公式
在这里插入图片描述

优点:精度高

ezplot('exp(x^2/2)');
hold on;
%四阶龙格-库塔方法
y = zeros(1,1000);
x = zeros(1,1000);
y(1)=1;
x(1)=0;
h=0.1;
for n=1:1000k1 = h*(x(n)*y(n));%h*f(x(n),y(n))k2 = h*((x(n)+h/2)*(y(n)+k1/2));%h*f(x(n)+h/2,y(n)+k1/2)k3 = h*((x(n)+h/2)*(y(n)+k2/2));%h*f(x(n)+h/2,y(n)+k2/2)k4 = h*((x(n)+h)*(y(n)+k3));%h*f(x(n)+h,y(n)+k3)y(n+1)=y(n)+(1/6)*(k1+2*k2+2*k3+k4);y(1)=1;x(n+1)=x(n)+h;
end
plot(x,y,'-');
xlim([0,6]);
ylim([0,100]);

在这里插入图片描述
细看是两条线
在这里插入图片描述

应用:小船渡河问题:

河宽100米,A(0,100) , O(0,0),水流速 V 1 V_1 V1,小船速度 V 2 V_2 V2,水流速度始终垂直OA并向下流冲(方向不变)
从O到A点,求轨迹路线图和时间
V 2 V_2 V2平行的相量 ( − x , 100 , − y ) (-x,100,-y) (x,100,y)

c o s θ = − x x 2 + ( 100 − y ) 2 cos_{\theta} = \dfrac{-x}{\sqrt{x^2}+(100-y)^2} cosθ=x2 +(100y)2x

s i n θ = 100 − y x 2 + ( 100 − y ) 2 sin_{\theta} = \dfrac{100-y}{\sqrt{x^2}+(100-y)^2} sinθ=x2 +(100y)2100y

{ d y d t = V 2 s i n θ d x d t = X 1 + V 2 c o s θ \begin{cases} \dfrac{dy}{dt}=V_2sin_{\theta}\\ \dfrac{dx}{dt}=X_1+V_2cos_{\theta} \end{cases} dtdy=V2sinθdtdx=X1+V2cosθ

d y d x = V 2 s i n θ V 1 + V 2 c o s θ = V 2 ( 100 − y ) x 2 + ( 100 − y ) 2 V 1 − V 2 x \dfrac{dy}{dx}=\dfrac{V_2sin_{\theta}}{V_1+V_2cos_{\theta}}=\dfrac{V_2(100-y)}{\sqrt{x^2+(100-y)^2}V_1-V_2x} dxdy=V1+V2cosθV2sinθ=x2+(100y)2 V1V2xV2(100y)

因为 x x x变量随着 y y y增大有增有减,而 y y y变量随着 x x x增大有递增,所以 y y y为自变量可以求解

变换后:
d x d y = V 1 + V 2 c o s θ V 2 s i n θ = x 2 + ( 100 − y ) 2 V 1 − V 2 x V 2 ( 100 − y ) \dfrac{dx}{dy}=\dfrac{V_1+V_2cos_{\theta}}{V_2sin_{\theta}}=\dfrac{\sqrt{x^2+(100-y)^2}V_1-V_2x}{V_2(100-y)} dydx=V2sinθV1+V2cosθ=V2(100y)x2+(100y)2 V1V2x

x 0 = y 0 = 0 x_0 = y_0 = 0 x0=y0=0

y = zeros(1,1000);
x = zeros(1,1000);
x(1)=0;
y(1)=0;
v1=30;
v2=200;
h=0.1;
%差分求解轨迹
for n=1:1000 %1000*h = 100x(n+1)=x(n)+h*((sqrt(x(n)^2+(100-y(n))^2)*v1-v2*x(n))/(v2*(100-y(n))));y(n+1)=y(n)+h;
endplot(x,y,'-');
hold on;
%四阶龙格-库塔方法求轨迹
y = zeros(1,1000);
x = zeros(1,1000);
x(1)=0;
y(1)=0;
for n=1:1000k1 = h*((sqrt(x(n)^2+(100-y(n))^2)*v1-v2*x(n))/(v2*(100-y(n))));%h*f(x(n),y(n))k2 = h*((sqrt((x(n)+h/2)^2+(100-(y(n)+k1/2))^2)*v1-v2*(x(n)+h/2))/(v2*(100-(y(n)+k1/2))));%h*f(x(n)+h/2,y(n)+k1/2)k3 = h*((sqrt((x(n)+h/2)^2+(100-(y(n)+k2/2))^2)*v1-v2*(x(n)+h/2))/(v2*(100-(y(n)+k2/2))));%h*f(x(n)+h/2,y(n)+k2/2)k4 = h*((sqrt((x(n)+h)^2+(100-(y(n)+k3))^2)*v1-v2*(x(n)+h))/(v2*(100-(y(n)+k3))));%h*f(x(n)+h,y(n)+k3)x(n+1)=x(n)+(1/6)*(k1+2*k2+2*k3+k4);y(n+1)=y(n)+h;
end
plot(x,y,'-');
%(sqrt(x(n)^2+(100-y(n))^2)
%差分求解时间
t = zeros(1,1000);
for n=1:1000t(n+1)=t(n)+h*(1/(v2*(100-y(n))/(sqrt(x(n)^2+(100-y(n))^2))));
endxlim([0,10]);
ylim([0,100]);

在这里插入图片描述
使用特解验证时间算的对不对
在这里插入图片描述在这里插入图片描述

进阶求二阶微分方程

例子:
{ y ′ ′ − 2 y ′ + y = 0 y 0 = 1 y 0 ′ = 2 \begin{cases} y''-2y'+y=0\\ y_0=1\\ y'_0=2 \end{cases} y′′2y+y=0y0=1y0=2
解析解为: y = e x + x e x y=e^x+xe^x y=ex+xex
差分求解:
y 1 = y , y 2 = y ′ y_1=y,y_2=y' y1=y,y2=y

{ y 1 ′ = y 2 y 2 ′ = 2 y 2 − y 1 y 1 0 = 1 y 2 0 = 2 \begin{cases} y_1'=y2\\ y_2'=2y_2-y_1\\ y_{1_0}=1\\ y_{2_0}=2 \end{cases} y1=y2y2=2y2y1y10=1y20=2

y1 = zeros(1,1000);
y2 = zeros(1,1000);
x = zeros(1,1000);
y1(1)=1;
y2(1)=2;
x(1)=0;
h=0.01;
xlim([0,6]);
ylim([0,100]);for n=1:1000y2(n+1)=y2(n)+h*(2*y2(n)-y1(n));y1(n+1)=y1(n)+h*y2(n);x(n+1)=x(n)+h;
end%plot(x,y1,'-');
hold on;
%四阶龙格-库塔方法
y1 = zeros(1,1000);
y2 = zeros(1,1000);
x = zeros(1,1000);
y1(1)=1;
y2(1)=2;
x(1)=0;
for n=1:1000k1 = h*(2*y2(n)-y1(n));%h*f(x(n),y(n))k2 = h*(2*(y2(n)+k1/2)-(y1(n)+h/2));%h*f(x(n)+h/2,y(n)+k1/2)k3 = h*(2*(y2(n)+k2/2)-(y1(n)+h/2));%h*f(x(n)+h/2,y(n)+k2/2)k4 = h*(2*(y2(n)+k3)-(y1(n)+h));%h*f(x(n)+h,y(n)+k3)y2(n+1)=y2(n)+(1/6)*(k1+2*k2+2*k3+k4);k1 = h*y2(n);%h*f(x(n),y1(n))k2 = h*(y2(n)+h/2);%h*f(x(n)+h/2,y(n)+k1/2)k3 = h*(y2(n)+h/2);%h*f(x(n)+h/2,y(n)+k2/2)k4 = h*(y2(n)+h);%h*f(x(n)+h,y(n)+k3)y1(n+1)=y1(n)+(1/6)*(k1+2*k2+2*k3+k4);x(n+1)=x(n)+h;
end
plot(x,y1,'-');%解析解ezplot('exp(x)+x*exp(x)');
xlim([0,6]);
ylim([0,100]);

在这里插入图片描述

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

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

相关文章

git多账号使用报错:You don‘t have permissions to push to “xxx/xxxx“ onGitHub. Would

git多账号使用报错:You don’t have permissions to push to “xxx/xxxx” onGitHub. Would 有的时候我们有两个甚至多个git账号(公司的git账号和自己的github),为了不混淆提交,我们需要在提交之前查看自己的git账号必…

9_1 Linux 网络管理

9_1 Linux 网络管理 文章目录 9_1 Linux 网络管理[toc]1. 构建yum仓库2. 配置主机名3. 配置网络参数之IP地址与子网掩码、网关地址3.1 修改网卡命名规则(eth0,eth1,eth2)3.2 配置地址3.2.1 nmcli3.2.2 配置文件修改IP地址、子网掩…

渗透测试之内核安全系列课程:Rootkit技术初探(五)

今天,我们来讲一下内核安全! 本文章仅提供学习,切勿将其用于不法手段! 目前,在渗透测试领域,主要分为了两个发展方向,分别为Web攻防领域和PWN(二进制安全)攻防领域。在…

人工智能对聊天机器人训练数据的“淘金热”可能会耗尽人类编写的文本

人工智能对聊天机器人训练数据的“淘金热”可能会耗尽人类编写的文本 像ChatGPT这样的人工智能系统可能很快就会耗尽让它们变得更聪明的东西——人们在网上写下和分享的数万亿字。 Epoch AI研究集团发布的一项新研究预计,科技公司将在大约十年之交——2026年至203…

Linux快速保存文件的快捷键:w和快速保存并退出ZZ

在Linux中,使用vi或Vim编辑器时,保存文件的快捷键主要是:w。以下是关于vi/Vim编辑器中保存操作的详细快捷键和说明: 1、保存文件 :w :这是最常用的保存命令。在命令模式下(按Esc键进入),输入:…

【Linux】ls命令

这个命令主要是用于显示指定工作目录下之内容(列出目前工作目录所含的文件及子目录)。 掌握几个重点的常使用的就可以: ls -l # 以长格式显示当前目录中的文件和目录 ls -a # 显示当前目录中的所有文件和目录&am…

Go使用https

一、服务端 1. 生成私钥和证书 安装OpenSSL windows安装OpenSSL生成CA证书创建证书 以上两个步骤,参考:Go http2 和 h2c 2. 代码 package mainimport ("log""net/http""time""golang.org/x/net/http2" )co…

TCP四次挥手全过程详解

TCP四次挥手全过程 有几点需要澄清: 1.首先,tcp四次挥手只有主动和被动方之分,没有客户端和服务端的概念 2.其次,发送报文段是tcp协议栈的行为,用户态调用close会陷入到内核态 3.再者,图中的情况前提是双…

【递归、搜索与回溯】穷举vs暴搜vs深搜vs回溯vs剪枝

穷举vs暴搜vs深搜vs回溯vs剪枝 1.全排列2.子集 点赞👍👍收藏🌟🌟关注💖💖 你的支持是对我最大的鼓励,我们一起努力吧!😃😃 管他什么深搜、回溯还是剪枝,画出决…

使用API有效率地管理Dynadot域名,创建文件夹管理域名

关于Dynadot Dynadot是通过ICANN认证的域名注册商,自2002年成立以来,服务于全球108个国家和地区的客户,为数以万计的客户提供简洁,优惠,安全的域名注册以及管理服务。 Dynadot平台操作教程索引(包括域名邮…

李廉洋:6.11黄金原油持续震荡上行,今日行情走势分析策略。

黄金消息面分析:上周黄金市场的走势受到了PCE通胀数据和美联储政策预期的显着影响。尽管市场对黄金的长期看涨情绪依然存在,但短期内金价的波动性预计将持续。4月份的PCE通胀数据显示价格压力有所降温,这一结果与分析师预期一致,但…

分享一些班组长日常管理工作的经验

在企业管理中,班组长作为基层的管理者和执行者,扮演着至关重要的角色。他们不仅要确保生产任务的顺利完成,还要负责团队建设、员工激励和沟通协调等工作。本文,深圳天行健精益管理咨询公司旨在分享一些班组长日常管理工作的经验&a…

探索C# 10.0的关键新特性

前言 随着.NET 6的发布,C# 10.0带来了许多创新特性和改进,旨在简化编码过程,增强开发者的生产力,并提供更现代、简洁的编程体验,可以说,这些新特性不仅增强了C#的表达能力,还提高了开发者的编程…

【c语言】文件操作,解开你的疑惑

文件操作 为什么使用文件什么是文件文件的分类文件名 二进制文件和文本文件文件的打开与关闭流与标准流流标准流 文件指针文件的打开与关闭 文件的顺序读写文件的随机读写文件读取结束的判定文件缓冲区 为什么使用文件 我们程序运行的数据是运行在内存中的,当成程序…

实用Python:文件与目录管理的17个技巧

今天我们要一起探索的是Python编程中的一个非常实用且基础的领域——文件与目录管理。无论是处理个人数据、自动化办公任务还是构建复杂的软件系统,这些技巧都将大大提升你的工作效率。准备好了吗?让我们一起动手吧! 1. 打开与读取文件 目标…

三生随记——梦魇之枕

第一章:搬入新家 在迷雾笼罩的小镇边缘,伫立着一座年代久远的木屋。李晴站在屋外,打量着这座看似平静却充满神秘感的居所。因为工作的原因,她不得不暂时搬离喧嚣的城市,来到这个陌生的地方。 木屋内部陈旧却别有一番风…

品牌渠道健康发展的关键与方法

一个品牌的渠道健康与否对其长期发展至关重要。品牌虽多,但并非所有产品都能成为品牌,创建品牌需大量精力,而让品牌长久健康发展则需多方面努力。 力维网络服务众多知名品牌,总结出一些渠道治理方法供品牌参考。首先,管…

【linux】(7)文本分析awk

awk 用于分析、过滤和生成报告。 基本用法 awk pattern {action} filename常用选项 -F:指定字段分隔符 awk -F, {print $1} filename例子:使用逗号作为字段分隔符,并打印第一列。 -v:定义变量 awk -v varvalue BEGIN {print va…

分享5款让大家电脑更好用的软件

​ 电脑是我们日常生活和工作中不可缺少的工具,今天给大家推荐了五款让电脑更好用的软件。 1.系统清理——CCleaner ​ CCleaner是一款系统优化和隐私保护工具,可以清理无用文件、浏览器缓存、回收站内容等,释放磁盘空间,提升系…

【linux网络(二)】网络基础之套接字编程

💓博主CSDN主页:杭电码农-NEO💓   ⏩专栏分类:Linux从入门到精通⏪   🚚代码仓库:NEO的学习日记🚚   🌹关注我🫵带你学更多操作系统知识   🔝🔝 Linux网络 1. 前言2. 端口号详…