最近在学计算机仿真,把一些简单的仿真方法分享一下,也帮自己记忆。
不是计算机专业,只是感兴趣选的选修课,代码写的不好,路过的大神轻喷。
- 基础欧拉公式
欧拉公式,一阶数值方法,可以说是入门级的了吧。简单介绍如下,其实就是根据当前点坐标求导数,乘以步长后加上当前的函数值即近似于下一个点的函数值,是一阶的,个人感觉在实用于导数变化很小的情况下导数变化大的情况下还是挺粗糙的,甚至有的时候还会跑飞。
对应的代码也挺简单的,如下:
#Eular Method
#欧拉公式,需要更改的地方用"####"标注
def eular_derivative(x0,y0):#返回x0,y0对应的导数值
return 0-y0**2####此处需要修改,为对应的导数值
def eular_method(start,end,step,y0):#返回给定条件下的欧拉公式的x,y值
y = [y0]
x = [start]
while x[-1]<end:
y.append(y[-1]+eular_derivative(x[-1],y[-1])*step)
x.append(x[-1]+step)
return x,y
eular_x,eular_y = eular_method(0,1,0.1,1)####此处需要修改x,y的初始值,步长,x的终值
- 改进的欧拉公式(improved Euler method)
欧拉公式其实不止前面所说的一种方式,出了上面所介绍的前向欧拉法外,还包括后向欧拉公式,中点欧拉公式等。但是这几种依旧是单步格式,即通过某一个点去推测另一个点的位置,本质上都类似,误差的可能性还是比较大的,他们也被成为显式欧拉公式,那么对应的,也就存在隐式欧拉公式。隐式欧拉公式其实就相当于将下一个点的值作为未知数代入方程中并求解,因此,其表达式如下:
代入求解,其精度必然较高,但是其计算量增大,当面对较复杂的方程时更是如此,求解方程的形式也给代码的编写带来了难度,因此,这种方式也有其限制条件。
到此,我们就可以很明显看到,隐式欧拉公式计算量小但精度有限,显式欧拉公式精度较高但是计算量大,那么,有没有什么方法可以综合其优缺点使得计算量与精度均在可接受范围内呢?
这就是改进的欧拉方法,也称为预测-矫正欧拉法。即先利用显式欧拉公式预测下一个值,再将下一个值导入隐式欧拉公式进行纠正,其误差减小,精度提高。
公式如下
通过公式,代码如下:
#improved Eular Method
#改进的欧拉公式,需要更改的地方用"####"标注
def derivative(x0,y0):#返回x0,y0对应的导数值
return 0-y0**2####此处需要修改,为对应的导数值
def improved_eular_method(start,end,step,y0):#返回给定条件下的改进的欧拉公式的x,y值
y = [y0]
x = [start]
while x[-1]<end:
Estimate_y = y[-1]+derivative(x[-1],y[-1])*step#预测值
y.append(y[-1]+(derivative(x[-1],y[-1])+Estimate_y)*step/2)#代入预测值求得精确值
x.append(x[-1]+step)
return x,y
improved_eular_x,improved_eular_y = improved_eular_method(0,1,0.1,1)####此处需要修改x,y的初始值,步长,x的终值
- 二者对比
那么,二者的精度等对比如何呢?下面举例进行对比.
以下面式子为例:
其精确值为
根据代码求解进行比较得:
好啦,哪个好不言而喻。
最后附上上面比较的代码,算是个代码总结:
#导数部分
def eular_derivative(x0,y0):
return 0-y0**2
#Eular Method
#欧拉公式
def eular_method(start,end,step,y0):
y = [y0]
x = [start]
while x[-1]<end:
y.append(y[-1]+eular_derivative(x[-1],y[-1])*step)
x.append(x[-1]+step)
return x,y
#improved Eular Method
#改进的欧拉公式
def improved_eular_method(start,end,step,y0):#返回给定条件下的改进的欧拉公式的x,y值
y = [y0]
x = [start]
while x[-1]<end:
Estimate_y = y[-1]+derivative(x[-1],y[-1])*step
y.append(y[-1]+(derivative(x[-1],y[-1])+derivative(x[-1],Estimate_y))*step/2)
x.append(x[-1]+step)
return x,y
#计算与输出
eular_x,eular_y = eular_method(0,1,0.1,1)#普通显式欧拉公式
improved_eular_x,improved_eular_y = improved_eular_method(0,1,0.1,1)#改进欧拉公式
len_num = len(eular_x)
print(" t\t 欧拉\t\t改进欧拉\t精确值")
for i in range(len_num):
Exact_solution = 1/(1+eular_x[i])#精确解
print(format(eular_x[i], '.2f')+"\t",end= '')
print(format(eular_y[i], '.6f')+" \t",end= '')
print(format(improved_eular_y[i], '.6f')+" \t",end ='')
print(format(Exact_solution, '.6f'))
第一次写文章,有用的话帮忙点个赞咯~