【计算方法】数值微积分

0.WARNINGS

本文章主要适用于计算方法代码的实现参考，由于本人是python究极小白，为了实验课速成了一些内容，因此会包含较多的暴力解法orz，如有代码错误、可优化的地方欢迎各位大佬指出，感激不尽。

此外，本人制作本文的目的，主要是懒的将代码保存在本地中，且博客方便个人的复习。（写这篇的时候已经知道下周要小测了，哇的一声哭出来！）

数值微积分主要包含的内容有：求积公式&代数精度；梯形公式；牛顿-科特斯公式；复化求积公式（Newton、Cotes、Gauss、Romberg、梯形等）；龙贝格公式；高斯型求积公式（切比雪夫、勒让德）；数值微分（向前、向后、中间）。

1.用辛普森公式求解定积分

import mathinp = input().split(' ')
a = float(inp[0])
b = float(inp[1])def fun(x):y=x*x*((math.e)**x)return ydef simpson(a,b):ret=1.0/6 * (b-a) * (fun(a) + 4*fun(a/2+b/2) + fun(b))return retprint("%.5f" % simpson(a,b))

2. 用科特斯公式计算定积分

import mathinp = input().split(' ')
a = float(inp[0])
b = float(inp[1])def fun(x):#x=0的时候分母没有意义！！！！！！！！！！！if x==0:y=1else:y=math.sin(x)/xreturn ydef cotes(a,b):temp=float((b-a)/4)ret=(b-a) * (fun(a)*(7/90) + fun(a+temp)*(16/45) + fun(a+temp*2)*(2/15) + fun(a+temp*3)*(16/45) +fun(b)*(7/90))return retresult = cotes(a,b)print("%.5f" % result)

3.比较梯形辛普森柯特斯三种不同的积分公式的求积结果

import mathinp = input().split(' ')
a = float(inp[0])
b = float(inp[1])def fun(x):y=x * math.sin(x) * math.e**xreturn ydef simpson(a,b):ret1=1.0/6 * (b-a) * (fun(a) + 4*fun(a/2+b/2) + fun(b))return ret1def cotes(a,b):temp=float((b-a)/4)ret2=(b-a) * (fun(a)*(7/90) + fun(a+temp)*(16/45) + fun(a+temp*2)*(2/15) + fun(a+temp*3)*(16/45) +fun(b)*(7/90))return ret2def ti(a,b):ret3=1.0/2 * (fun(a) + fun(b)) * (b-a)return ret3myti=ti(a,b)
mysimpson=simpson(a,b)
mycotes=cotes(a,b)
print("%.5f %.5f %.5f" % (myti,mysimpson ,mycotes))

4. 用复化梯形公式计算式子的近似值

import mathinp=input().split(' ')
a=float(inp[0])
b=float(inp[1])
n=int(inp[2])def fun(x):y=math.log(x)return ydef ti(a,b,n):h=(b-a)/nT=0for i in range(1,n):T += fun(a+i*h)T = h/2 * (fun(a)+2*T+fun(b))return Tret=ti(a,b,n)
print("%.5f" % ret)

5. 用复化辛卜生公式计算积分的近似值

import mathinp=input().split(' ')
a=float(inp[0])
b=float(inp[1])
n=int(inp[2])def fun(fai):y=math.sqrt(4-math.sin(fai)*math.sin(fai))return ydef simpson(a,b,n):h=(b-a)/nS=0for i in range(1,n):S += 2*fun(a+i*h)j=a+h/2while j<b:S += 4*fun(j)j += hS = h/6 * (fun(a)+S+fun(b))return Sret=simpson(a,b,n)
print("%.5f" % ret)

6.变步长梯形积分求解定积分

import mathinp=input().split(' ')
a=float(inp[0])
b=float(inp[1])
e=float(inp[2])def fun(x):if x==0:return 1else:y=math.sin(x)/xreturn yh=b-a
T1 =h/2*(fun(b)+fun(a))while 1:s=0x=a+h/2while(x<b):s+=fun(x)x+=hT2=T1/2+h/2*sif(math.fabs(T2-T1)<=e):breakelse:T1=T2h/=2ret=T2
print("%.6f" % ret)

7. 用龙贝格积分法计算积分

import math
import numpy as npinp=input().split(' ')
a=float(inp[0])
b=float(inp[1])def fun(x):y=4/(1+x**2)return ydef romberg(a,b,n):cnt=0while n>1:n/=2cnt+=1sum=cnt+1 + 3array=np.zeros(sum)# T1 is array[0]h=b-aarray[0]=h/2 * (fun(a)+fun(b))# print(array[0])for i in range(1,sum):s=0x=a+h/2while x<b:s+=fun(x)x+=hT2=array[i-1]/2+h/2*sh/=2array[i]=T2# print(array[i])# Sncount=sum-1# print(count)while count>0:array[count]=array[count]*4/3-array[count-1]*1/3# print(array[count])count -= 1# Cncount=sum-1while count>1:array[count]=array[count]*16/15-array[count-1]*1/15# print(array[count])count -= 1# Rncount=sum-1while count>2:array[count] = array[count] * 64 / 63 - array[count - 1] * 1 / 63# print(array[count])count -= 1# output R4print("%.2f" % array[sum-1])
"""# from R1 ~ Rnfor i in range(3,sum):print(array[i])
"""romberg(a,b,4)

8.用n=2高斯-勒让德求积公式计算下列积分

import mathinp = input().split(' ')
a = float(inp[0])
b = float(inp[1])def fun(x):y=math.sqrt(x+1.5)return ydef gauss(a,b):x0=-math.sqrt(15)/5x1=0x2=math.sqrt(15)/5myx0=(b-a)/2 * x0 + (a+b)/2myx1 = (b - a) / 2 * x1 + (a + b) / 2myx2 = (b - a) / 2 * x2 + (a + b) / 2I=5/9 * fun(myx0)+ 8/9*fun(myx1 )+5/9*fun(myx2)I*=(b-a)/2return Iret=gauss(a,b)
print("%.5f" % ret)

9. 比较4种高精度的求积方法

import math
import numpy as npinp = input().split(' ')
a = float(inp[0])
b = float(inp[1])def fun(x):if x==0:return 0else:return 1/xdef romberg(a,b,n):cnt=0while n>1:n/=2cnt+=1sum=cnt+1 + 3array=np.zeros(sum)# T1 is array[0]h=b-aarray[0]=h/2 * (fun(a)+fun(b))# print(array[0])for i in range(1,sum):s=0x=a+h/2while x<b:s+=fun(x)x+=hT2=array[i-1]/2+h/2*sh/=2array[i]=T2# print(array[i])# Sncount=sum-1# print(count)while count>0:array[count]=array[count]*4/3-array[count-1]*1/3# print(array[count])count -= 1# Cncount=sum-1while count>1:array[count]=array[count]*16/15-array[count-1]*1/15# print(array[count])count -= 1# Rncount=sum-1while count>2:array[count] = array[count] * 64 / 63 - array[count - 1] * 1 / 63# print(array[count])count -= 1# return R4return array[sum-1]def change(a,b,x):x=(b-a)/2 * x + (a+b)/2return xdef twolegendre(a,b):I=0x0=-1/math.sqrt(3)x1=-x0x0=change(a,b,x0)x1=change(a,b,x1)I=fun(x0)+fun(x1)I*=(b-a)/2return Idef threelegendre(a,b):I=0x0=-math.sqrt(15)/5x1=0x2=-x0x0=change(a,b,x0)x1=change(a,b,x1)x2 = change(a, b, x2)I=5/9*fun(x0)+8/9*fun(x1)+5/9*fun(x2)I*=(b-a)/2return Idef fivelegendre(a,b):I=0x0=-0.9061798x1=-0.5384693x2=0x3=-x1x4=-x0x0 = change(a, b, x0)x1 = change(a, b, x1)x2 = change(a, b, x2)x3 = change(a, b, x3)x4 = change(a, b, x4)A0=0.2369269A1=0.4786287A2=0.5688889I=A0*(fun(x0)+fun(x4))+A1*(fun(x1)+fun(x3))+A2*fun(x2)I *= (b - a) / 2return Idef complexgauss(a,b):I=0h=(b-a)/4for i in range(0,4):newa=a+i*hnewb=a+(i+1)*hI+=twolegendre(newa,newb)return Iret1=romberg(a,b,4)
ret2=threelegendre(a,b)
ret3=fivelegendre(a,b)
ret4=complexgauss(a,b)
print("%.5f" % ret1)
print("%.5f" % ret2)
print("%.5f" % ret3)
print("%.5f" % ret4)

10.用高斯切比雪夫求积公式求积分

import math
"""
inp = input().split(' ')
a = float(inp[0])
b = float(inp[1])
"""
def fun(x):y=math.e**xreturn ydef chebyshev(n):A=math.pi/(n+1)I=0for i in range(0,n+1):x=math.cos((2*i+1)/(2*n+2)*math.pi)I+=fun(x)I*=Areturn In=int(input())
# x0,x1,...xn个高斯点
ret=chebyshev(n)
print("%.6f" % ret)

11.用向前差商公式、向后差商公式和中心差商公式计算定在sinx/x点的导数

import math
"""
inp = input().split(' ')
a = float(inp[0])
b = float(inp[1])
"""
def fun(x):if x==0:return 1else:y = math.sin(x) / xreturn ydef forward(x):h=1temp=0while 1:f=1/h * (fun(x+h)-fun(x))if math.fabs(f-temp)<0.001:breakelse:temp=fh/=2return temp # 输出前点的导数值def backward(x):h=1temp=0while 1:f=1/h * (fun(x)-fun(x-h))if math.fabs(f-temp)<0.001:breakelse:temp=fh/=2return tempdef mid(x):h = 1temp = 0while 1:f=1/2/h * (fun(x+h)-fun(x-h))if math.fabs(f-temp)<0.001:breakelse:temp=fh/=2return tempmyx=float(input())
print("%.3f %.3f %.3f" % (forward(myx),backward(myx),mid(myx)))

12.利用数值求导的三点公式计算人口增长率

import math
import numpy as npyear=[float(e) for e in input().split()]population=[float(t) for t in input().split()]length=len(year)
# print(length)h=year[1]-year[0]
# print(h)
rate=np.zeros(length)A=1/(2*h)
for i in range(0,length):if i==0:rate[i]=A*(-3*population[0]+4*population[1]-population[2])elif i==length-1:rate[i]=A*(population[length-3]-4*population[length-2]+3*population[length-1])else:rate[i]=A*(population[i+1]-population[i-1])print("%.3f" % rate[i])

【计算方法】数值微积分相关推荐

工程数学（数值分析）第六讲：数值微积分
文章目录第六讲:数值微积分二点.三点数值微积分求导数定步长梯形求积公式变步长梯形求积公式定步长.变步长辛普森求积公式第六讲:数值微积分二点.三点数值微积分求导数定步长梯形求积公式变步 ...
【台大郭彦甫】Matlab入门教程超详细学习笔记七：数值微积分（附PPT链接）
数值微积分前言一.多项式微积分 1. 多项式计算 2. 多项式微分 3. 多项式积分二.数值微积分 1. 数值微分法 2. 高阶微分法 3. 数值积分法三.回顾Function Handles ...
matlab武汉理工大学数值分析线性函数拟合实验_11数值分析第七章数值微积分龙贝格积分大学数学云课堂...
数值分析基础入门所有知识点相关内容,请按照前面数字依次学习 1第二章数值分析的基本概念之有效数字与秦九韶算法 1数值分析第二章基本概念作业答疑大学数学云课堂 2数值分析第三章线性代数方程组的直接法大学 ...
matlab的数值计算方法,数值计算方法中的一些常用算法的Matlab源码
数值计算方法中的一些常用算法的Matlab源码,这些程序都是原创,传上来仅供大家参考,不足之处请大家指正,切勿做其它用途-- 说明:这些程序都是脚本函数,不可直接运行,需要创建函数m文件,保存时文件名 ...
MATLAB设置x为0到10所有数,MATLAB教学_10数值微积分
本文学习视频地址:https://www.bilibili.com/video/av68228488?p=10 课堂PPT以及本人学习代码已上传. 本文学习内容: 多项式的微分和积分数值的微分和积分 ...
数值微积分与方程求解
1.数值微分与数值积分数值微分 MATLAB提供了求向前差分的函数diff,其调用格式有三种: dx=diff(x):计算向量x的一阶向前差分,dx(i)=x(i+1)-x(i),i=1,2,-,n ...
matlab trapz二重积分函数_matlab数值微积分
1.polyval() %多项式构造函数,参数为系数vector,自变量vector f=[9,-5,3,7]; x=-2:0.01:5; %x的范围为-2到5 y=polyval(f,x ...
MATLAB（九）数值微积分
文章目录前言用MATLAB表示多项式多项式的值:polyval() 多项式微分:polyder() Exercise练习多项式的集成:polyint() 偏差:diff() Exercise练 ...
数学建模之matlab软件学习06——专题六数值微积分与方程求解
#本文章仅用于记录本人学习过程,当作笔记来用,如有侵权请及时告知,谢谢! 6.1 数值微分与数值积分数值微分 6.2 6.3 线性方程组应用举例: 平面桁架结构受力分析问题小行星运行轨道计算问题: ...

【计算方法】数值微积分

【计算方法】数值微积分相关推荐

最新文章

热门文章