Python数值求解微分方程方法(欧拉法,隐式欧拉)
作者:Chandler_river 发布时间:2023-06-29 10:45:29
标签:Python,数值,求解
不说什么,先上代码
这里先求解形如的微分方程
1.欧拉法
def eluer(rangee,h,fun,x0,y0):
step = int(rangee/h)
x = [x0] + [h * i for i in range(step)]
u = [y0] + [0 for i in range(step)]
for i in range(step):
u[i+1] = u[i] + h * fun(x[i],u[i])
plt.plot(x,u,label = "eluer")
return u
2.隐式欧拉法
def implicit_euler(rangee,h,fun,x0,y0):
step = int(rangee/h)
x = [x0] + [h * i for i in range(step)]
u = [y0] + [0 for i in range(step)]
v = ["null"] + [0 for i in range(step)]
for i in range(step):
v[i+1] = u[i] + h * fun(x[i],u[i])
u[i+1] = u[i] + h/2 * (fun(x[i],u[i]) + fun(x[i],v[i+1]))
plt.plot(x,u,label = "implicit eluer")
return u
3.三阶runge-kutta法
def order_3_runge_kutta(rangee,h,fun,x0,y0):
step = int(rangee/h)
k1,k2,k3 = [[0 for i in range(step)] for i in range(3)]
x = [x0] + [h * i for i in range(step)]
y = [y0] + [0 for i in range(step)]
for i in range(step):
k1[i] = fun(x[i],y[i])
k2[i] = fun(x[i]+0.5*h,y[i]+0.5*h*k1[i])
k3[i] = fun(x[i]+0.5*h,y[i]+2*h*k2[i]-h*k1[i])
y[i+1] = y[i] + 1/6 * h * (k1[i]+4*k2[i]+k3[i])
plt.plot(x,y,label = "order_3_runge_kutta")
return y
4.四阶runge-kutta法
def order_4_runge_kutta(rangee,h,fun,x0,y0):
step = int(rangee/h)
k1,k2,k3,k4 = [[0 for i in range(step)] for i in range(4)]
x = [x0] + [h * i for i in range(step)]
y = [y0] + [0 for i in range(step)]
for i in range(step):
k1[i] = fun(x[i],y[i])
k2[i] = fun(x[i]+0.5*h,y[i]+0.5*h*k1[i])
k3[i] = fun(x[i]+0.5*h,y[i]+0.5*h*k2[i])
k4[i] = fun(x[i]+h,y[i]+h*k3[i])
y[i+1] = y[i] + 1/6 * h * (k1[i]+2*k2[i]+2*k3[i]+k4[i])
plt.plot(x,y,label = "order_4_runge_kutta")
return y
5.上图
当然,想要成功操作,得加上这个
rangee = 1
fun = lambda x,y:y-2*x/y
implicit_euler(rangee,0.0001,fun,0,1)
order_4_runge_kutta(rangee,0.0001,fun,0,1)
order_3_runge_kutta(rangee,0.0001,fun,0,1)
eluer(rangee,0.0001,fun,0,1)
plt.legend()
plt.show()
来源:https://blog.csdn.net/Chandler_river/article/details/124295877
0
投稿
猜你喜欢
- 一、前言我们在使用pycharm写代码时可能出现过下面这种情况,不小心点到了ignore ....:这样会导致整个代码都没有错误提示了,类似
- 当然,页面最好不要刷新,但是,拷贝一下浏览器的链接,又希望是下次能定位到你播发的那个视频。方法很简单,改变一下 url 的 fragment
- 说起计算机中的时间,还有一些比较有意思的事,比如我们经常听到的Unix时间戳,UTC时间,格林威治时间等,从表示上来讲他们基本属于同一个东西
- 对数据库的管理常规就是进行预防性的维护,以及修复那些出现问题的内容。进行检查和修复通常具有四个主要的任务:1. 对表进行优化2. 对表进行分
- 说实话,第一次写这么长的Python代码,期间遇到了很多问题,但是,最终还是完成了,花了我一天半的时间。该程序实现了用户的增,删,改,查,主
- 目录引言环境代码实现准备工作获取并存储好友头像信息生成微信好友墙后记引言前段时间,微信朋友圈开始出现了一种晒照片新形式,微信好友墙,即在一张
- 在用户体验量化方法研究(一)中,我们以用户体验层次模型为基础,提出了三种量化用户体验的方法:以任务为中心、以行为为中心以及以体验为中心的量化
- 我们先来实现一个简单的例子,hello world。似乎每种语言教程的第一节都会讲这个,我们也不例外。首先我们先创建一个项目目录,目录可自己
- JS无法进行精确计算的bug在做CRM,二代审核需求审核详情页面时。需要按比例(后端传类似0.8的小数)把用户输入的数字显示在不同的地方。
- 匹配开头结尾字符功能^匹配字符串开头$匹配字符串结尾示例1:$需求:匹配163.com的邮箱地址email_list = [ "
- Windows Registry Editor Version 5.00 [HKEY_LOCAL_MACHINE\SOFTWARE\Micr
- 一、介绍对于Visual Studio Code开发工具,有一款优秀的GoLang插件,它的主页为:https://github.com/m
- Python django中我们经常用的response有django中的 JsonResponse, HttpResponse,还有DRF
- 本文实例讲述了Python简单读写Xls格式文档的方法。分享给大家供大家参考,具体如下:1. 模块安装使用pip install命令安装,即
- 分组取TOP数据是T-SQL中的常用查询, 如学生信息管理系统中取出每个学科前3名的学生。这种查询在SQL Server 2005之前,写起
- 本文实例讲述了python获取本地计算机名字的方法。分享给大家供大家参考。具体如下:import sys, sockethostname =
- 分析使用CrawlSpider结合LinkExtractor和Rule爬取网页信息LinkExtractor用于定义链接提取规则,一般使用a
- tkinter介绍tkinter是python自带的GUI库,是对图形库TK的封装tkinter是一个跨平台的GUI库,开发的程序可以在wi
- 最近没有项目做,闲来无事写了一个小demo,特此分享到脚本之家平台,供大家参考下,本文写的不好还请各位大侠见谅!功能及方法逻辑都注释在代码中
- print() 方法用于打印当前窗口的内容,支持部分或者整个网页打印。调用 print() 方法所引发的行为就像用户单击浏览器的打印按钮。通