一阶微分方程的求解

算法的数学描述图解

wns威斯尼斯人 1

以下只是例题介绍,编程方式情参照例题,因为markdown用的不好,现在我都不知道怎样输入数学公式。所以私底下还得记笔记。今天偶然读了一篇文章种,文章中说,如今的大学者越发少见。其中一个原因是现金我们重理轻文。的确如此啊,理学学科应用的太广泛,但是我还是学的是皮毛,虽然很懊恼。但是还得耐着性子一步步的慢慢来。大神不是一日练成的哦!

本篇介绍一下一阶微分方程的求解方法,以及伯努利方程的特殊求解方法。这个应该是上学时高数课中的内容,现在用到了,温习一下。
顺便感叹一下,时间过得真快。

wns威斯尼斯人 2

python求解水仙花数的方法,python求解水仙花

本文实例讲述了python求解水仙花数的方法。分享给大家供大家参考。具体如下:

一个N位的十进制正整数,如果它的每个位上的数字的N次方的和等于这个数本身,则称其为花朵数。

#!/usr/bin/python
def get_flower(n, ofile):
 D_pow=[pow(i,n) for i in range(0,10)]
 V_min=1*pow(10,n-1)
 V_max=sum((9*pow(10,x) for x in range(0,n)))
 T_count=0
 print D_pow, V_max, V_min
 nums=[1]+[0]*(n-1)
 print 'Start:', nums
 idx=n-1
 tmp_l=[0]*10
 while True:
  nums[idx]+=1
  if nums[idx]<10:
   j=idx+1
   while j<n:
    nums[j]=nums[idx] # reset 
    j+=1
   v=sum((D_pow[x] for x in nums))
   if v<=V_max and v>=V_min:
    T_count+=1
    #test if is flower
    #print 'do test:', ''.join(map(str,nums))
    k=0
    while k<10:
     tmp_l[k]=0
     k+=1
    N=0
    for k in nums:
     tmp_l[k]+=1
     N+=1
    while N>0:
     p=v%10
     if tmp_l[p]>0:
      tmp_l[p]-=1
      N-=1
     else:
      break
     v/=10
    if N==0:
     print >>ofile, 'hit', sum((D_pow[x] for x in nums))
   idx=n-1
  elif idx==0:
   print 'done'
   break
  else:
   idx-=1
 print 't_count', T_count
if __name__ == '__main__':
 with file('./f.txt', 'wb') as o:
  get_flower(21, o)
  #get_flower(3, o)

希望本文所述对大家的Python程序设计有所帮助。

本文实例讲述了python求解水仙花数的方法。分享给大家供大家参考。具体如下:
一个N位的十进…

实例

用Euler算法求解初值问题
\[ \frac{dy}{dx}=y+\frac{2x}{y^2}\]
初始条件\,自变量的取值范围\(x \in [0, 2]\)

例1:一阶微分方程

y=dsolve(‘Dy=1+y^2′,’y(0)=1′,’x’)
y =tan(pi/4 + x)

1. 定义

wns威斯尼斯人 3

形如上式的方程称为一阶线性微分方程,
并且当Q(x)恒为零时称为齐次线性方程,
Q(x)不恒为零时称为非齐次线性方程.

 

一、代数多项式法:

算法Python3代码求解

# 导入包import numpy as npimport matplotlib.pyplot as plt# 定义求解函数 y_dot = y + 2*x/def fx:    return y + 2*x/# 算法定义def ode_euler(f, y0, tf, h):    """    Solve and ODE using Euler method.    Solve the ODE y_dot = f    Parameters    ------------    :param f: function            Function describing the ODE    :param y0: array_like            Initial conditions.    :param tf: float            Final time.    :param h: float            Time step    :return:    y : array_like        Solution to the ODE.    t : array_like        Time vector.    """    y0 = np.array    ts = np.arange(0, tf + h, h)    y = np.empty((ts.size, y0.size))    y[0, :] = y0    for t, i in zip(ts[1:], range(ts.size - 1)):        y[i + 1, :] = y[i, :] + h * f(y[i, :], t)    return y, ts# 实例应用案例def newton_cooling_example():    print('Solving Newton Cooling ODE...')    y, ts = ode_euler(fx, 1, 2, 0.01)    print('Done.')    plt.figure()    plt.plot    plt.xlabel('time [s]')    plt.title('Solution to the Newton cooling equation')    plt.show()

例2:常系数的二阶微分方程

dsolve(‘D2y-2Dy-3y=0′,’y(0)=1,Dy(0)=0′,’x’)
ans =
(3exp(-x))/4 + exp(3x)/4

2. 通解

 1 tic;
 2 clear
 3 clc
 4 % N=input('please key in the value of ''N''');
 5 N=10;
 6 M=100;
 7 h=1/M;
 8 X=0:h:1;
 9 accurate_fun=inline('x.^2 - (2*exp(x))/(exp(1) + 1) - (2*exp(-x)*exp(1))/(exp(1) + 1) + 2');
10  f=inline('x.^2-x');
11  phi=inline('x.*(1-x).*x.^(i-1)','i','x');
12  diff_phi=inline('i*x.^(i-1)-(i+1)*x.^i','i','x');
13  for j=1:N
14      for i=1:N
15 A(i,j)=quad(@(x)phi(i,x).*phi(j,x)+diff_phi(i,x).*diff_phi(j,x),0,1);
16      end
17      b(j,1)=quad(@(x) phi(j,x).*f(x),0,1);
18  end
19  C=A\b;
20 syms x;
21  Un=0;
22 for i=1:N
23 Un=Un+C(i)*phi(i,x);
24 end
25 Un=Un+x;
26  numerical= double(vpa(subs(Un,'x',X)));
27  accurate=accurate_fun(X);
28  subplot(1,2,1)
29  plot(X,numerical,'r -',X,accurate,'b >');
30   title('numerical VS accurate');
31  legend('numerical','accurate');
32  grid on;
33   subplot(1,2,2);
34   plot(X,numerical-accurate,'g');
35   title('error');
36   grid on;
37  toc;

代码中的部分函数理解

例3:非常系数的二阶微分方程

dsolve(‘D2x-(1-x^2)*x+x=0′,’x(0)=3,Dx(0)=0’)
ans =
[ empty sym ]

无解咋办,欲知详解请看下文哦!

2.1 齐次线性方程的通解

对于齐次线性方程:

wns威斯尼斯人 4

可以推出:

wns威斯尼斯人 5
通解为:

wns威斯尼斯人 6奥门威胁斯wns威斯尼斯人 , 

 

numpy.array

numpy.array(object, dtype=None, copy=True, order='K', subok=False, ndmin=0)
参考numpy.array
output:创建一个array,返回类型为ndarray
实例

np.array([1, 2, 3.0]) # array([1., 2., 3.])np.array([[1, 2], [3, 4]]) # array([[1, 2], [3, 4]])np.array([1, 2, 3], dtype=complex) # array([1.+0.j, 2.+0.j, 3.+0.j])

例4:非线性微分方程

x=dsolve(‘(Dx)2+x2=1′,’x(0)=0′)

x =

cosh((pii)/2 + ti)
cosh((pii)/2 – ti)

2.2  非齐次线性方程的通解

对于非齐次线性方程:

wns威斯尼斯人 7设通解为:

wns威斯尼斯人 8

带入非齐次线性方程:

wns威斯尼斯人 9积分得:

wns威斯尼斯人 10其中C为常数。 

于是非齐次线性通解是:

wns威斯尼斯人 11

由此可以看出,齐次线性方程的通解是非齐次线性方程的一个特解。

 

wns威斯尼斯人 12

numpy.arange

参考numpy.arange
numpy.arange([start, ]stop, [step, ]dtype=None)
作用:在给定间隔内返回均匀间隔的值。
值在半开区间[start,
stop)内生成(换句话说,包括开始但不包括终止)。返回的是ndarray而不是列表。
np.arange()函数返回一个有终点和起点的固定步长的排列,如[1,2,3,4,5],起点是1,终点是5,步长为1。
参数个数情况: np.arange()函数分为一个参数,两个参数,三个参数三种情况

1. 一个参数时,参数值为终点,起点取默认值0,步长取默认值1。 2. 两个参数时,第一个参数为起点,第二个参数为终点,步长取默认值1。 3. 三个参数时,第一个参数为起点,第二个参数为终点,第三个参数为步长。其中步长支持小数。

案例

np.arange # array([3, 4, 5, 6])np.arange # array

例5:微分方程组

[x,y]=dsolve(‘Dx=3x+4y’,’Dy=-4x+3y’,’x(0)=0,y(0)=1′)
x =
sin(4t)exp(3t)
y =
cos(4
t)exp(3t)

3.  伯努利方程

wns威斯尼斯人 13

形如上式的方程叫做伯努利方程。

将方程线性化得:

wns威斯尼斯人 14

例子:

求下列方程的通解

wns威斯尼斯人 15两边除以y的平方得:

wns威斯尼斯人 16将第一项中y的负平方移入微分内得:

wns威斯尼斯人 17由非齐次线性方程的通解可知:

wns威斯尼斯人 18即原方程的通解为:

wns威斯尼斯人 19

 

numpy.ma.size

numpy.ma.size(obj, axis=None)
参考
案例

a = np.array([[1,2,3],[4,5,6]])np.size # 6np.size # 3np.size # 2

没有解析解可以采用数值解的方法

[x,y]=ode23(‘weif’,[0,1],1)
plot(x,y,’r’);
hold on
ezplot(‘x+exp(-x)’,[0,1])
x =
0
0.1000
0.2000
0.3000
0.4000
0.5000
0.6000
0.7000
0.8000
0.9000
1.0000
y =
1.0000
1.0048
1.0187
1.0408
1.0703
1.1065
1.1488
1.1966
1.2493
1.3066
1.3679

二、三角函数法:

numpy.empty

参考
numpy.empty(shape, dtype=float, order='C')
shape : int or tuple of int Shape of the empty array, e.g., or 2.
out : ndarray
案例

np.empty# 结果array([[ -9.74499359e+001, 6.69583040e-309],       [ 2.13182611e-314, 3.06959433e-309]]) #randomnp.empty([2, 2], dtype=int)# 结果array([[-1073741821, -1067949133],       [ 496041986, 19249760]]) #random

van der pol方程(就是例3,数值法求解)

[t,y]=ode23(‘vdpol’,[0,20],[3,0]);
y1=y(:,1);
y2=y(:,2);
plot(t,y1,’b’,t,y2,’r–‘)
legend(‘y1′,’y2’)

wns威斯尼斯人 20

van der pol

看看你能不能提出问题呢?所以我不吧我的疑惑写出来。嘿嘿

 1 tic;
 2 clear
 3 clc
 4 % N=input('please key in the value of ''N''');
 5 N=10;
 6 M=100;
 7 h=1/M;
 8 X=0:h:1;
 9 accurate_fun=inline('x.^2 - (2*exp(x))/(exp(1) + 1) - (2*exp(-x)*exp(1))/(exp(1) + 1) + 2');
10  f=inline('x.^2-x');
11  phi=inline('sin(i*pi*x)','i','x');
12  diff_phi=inline('i*pi*cos(i*pi*x)','i','x');
13  for j=1:N
14      for i=1:N
15 A(i,j)=quad(@(x)phi(i,x).*phi(j,x)+diff_phi(i,x).*diff_phi(j,x),0,1);
16      end
17      b(j,1)=quad(@(x) phi(j,x).*f(x),0,1);
18  end
19  C=A\b;
20  syms x;
21  Wn=0;
22 for i=1:N
23 Wn=Wn+C(i)*phi(i,x);
24 end
25 Un=Wn+x;
26 numerical=vpa(subs(Un,'x',X));
27 accurate=accurate_fun(X);
28 subplot(1,2,1)
29 plot(X,numerical,'r -',X,accurate,'g ^');
30 title('Ritz Galerkin method');
31 legend('numerical solution','accurate solution');
32 grid on;
33 subplot(1,2,2)
34 plot(X,numerical-accurate,'b');
35 title('error');
36 grid on;
37 toc;

 

 wns威斯尼斯人 21

 

相关文章