用有限差分和牛顿法解非线性微分方程（边值问题）-python

最近师兄让我解一个微分方程，我随口就答应下来了。结果仔细研究以后发现是个大坑。方程式一个复杂的非线性方程，而且是边值问题（知道两个端点的值）。微分方程的初值问题（知道一个端点的值和导数）相对简单，因为可以降阶，但边值问题不能降阶，相对麻烦一些，先上方程：
$(1-{\mu }^2)T^{\prime \prime }(\mu )+v(v-1)\sqrt{T(\mu )}-(1-{\mu }^2{)}^{2}T(\mu )T^{\prime \prime }(\mu )M^{2}T(\mu {)}^{-1-2/v}=0$
边界条件是
$T(1)=0和T(0)=1$
这个方程我没解出来，知道好的算法的小伙伴麻烦在下面留言.
如果只有前面两项，会容易一点
$(1-{\mu }^2)T^{\prime \prime }(\mu )+v(v-1)\sqrt{T(\mu )}=0$
v是参数，取1的时候有解析解$T=1-\mu$
边值问题一般用有限差分法和打靶法，这里我用的是有限差分。
在$\mu$区间[0,1]上分成n等分，步长h，T在i点的二阶导数可以表示成
$T_i^{\prime \prime }=(T_{i+1}-2T_i+T_{i-1})/h^2$
T在每个点上的值作为未知数，这样每个点上都有一个未知数，一共有n个点，端点已知，就是n-2个未知数。把上面这个式子带到待求的微分方程里面，每个点上列一个方程，有n-2个方程。如果是线性方程，可以直接用追赶法，高斯-赛的尔法等来解，由于这个方程是非线性的，这里我用牛顿法来解。关于牛顿法见我上一篇文章：python 实现（拟）牛顿法解非线性方程组
https://blog.csdn.net/zry1318/article/details/86033518

import math
import numpy as np  
from numpy import * 
import matplotlib.pyplot as plt

def Fun(x,num):                       #方程组在这里
    i  =  num
    f  =  np.zeros((i),dtype = float)
    ddT0 = (x[1]-2*x[0]+T0)/h/h#  待求函数二阶导数，先列靠近左边界点的方程
    muj10 = 1.-mu[1]*mu[1] # 1-mu*mu
    f[0]  = muj10*ddT0+vi*(vi-1.)*x[0]**0.5

    ddTn = (T1-2*x[i-1]+x[i-2])/h/h#  待求函数二阶导数，靠近右边界点的方程
    muj1n = 1.-mu[i]*mu[i] # 1-mu*mu
    f[i-1]  = muj1n*ddTn+vi*(vi-1.)*x[i-1]**0.5

    for j in range(1, i-1):#边界点知道，j从1开始
        ddT = (x[j+1]-2*x[j]+x[j-1])/h/h#  待求函数T(mu)的二阶导数
        muj1 = 1.-mu[j]*mu[j] # 1-mu*mu

        f[j]  =  muj1*ddT+vi*(vi-1.)*x[j]**0.5
    return f

def dfun(x,num):                         #计算雅可比矩阵的逆矩阵
    df  =  np.zeros((num,num),dtype = float)
    dx  =  0.000001                           # 
    x1  =  np.copy(x)
    for i in range(0,num):              # 求导数，i是列，j是行
        for j in range(0,num):
            x1  =  np.copy(x)
            x1[j]  =  x1[j]+dx           #x+dx
            df[i,j]  =  (Fun(x1,num)[i]-Fun(x,num)[i])/dx   #f(x+dx)-f（x）/dx
    df_1  =  np.linalg.inv(df)                              #计算逆矩阵
    return df_1

def Newton(x,num):
               
    x1  =  np.copy(x)
    i  =  0
    delta  =  np.copy(x)
#    dfun0 = dfun(x,num)          #也可以使用简化牛顿法
    while(np.max(abs(delta)) > 1.e-6 and i < imax):  #控制循环次数
        print(i)
        x1  =  x-dot(dfun(x,num),Fun(x,num))  #公式
        mask=x1<0
        x1[mask]=x_stand[mask]
        delta  =  x1-x                     #比较x的变化
        x  =  x1
        i  =  i+1
        print(x)
    if i == imax:
        print('error: i>imax, no solution')
    print('i = ',i)
    return x


mu = np.arange(0.0, 1.01,  0.05)   #自变量
h = mu[1]-mu[0]#步长，线性
MM = 0.0#   M^2
vi =  1.
T0 = 1.#边界条件,
T1 = 0.
imax = 10000#最大迭代次数
yy = 1.-mu    #精确解
x_stand = yy[1:-1]
num  = mu.size-2                      # 方程未知数的个数
x  =  np.ones((num),dtype = float)/2.#初始值
a  =  Newton(x,num)

print('imax = ',imax)
print('note: if i  =  imax, the code did not find solution')
print('final result = ',a)
plt.figure(0) # 创建图表1    

#print(yy)  
a2  =  np.zeros((mu.size),dtype = float)
a2[0] = 1
a2[1:-1] = a

plt.plot(mu, a2,'r' ) #绘制数据散点图，用plt.sc  atter(x,y)
print(a2)
print(yy)
plt.plot(mu, yy,'b' ) #绘制数据散点图，用plt.sc  atter(x,y)

plt.xlabel('$\mu$')     #设置X轴的文字   
plt.ylabel('T')  
#plt.savefig('Figure_T.png') 
plt.show() #显示图片

v=1的时候，数值解与解析解重合

v=2的时候，结果如下