SciPy笔记（有初始条件的常微分方程）
Posted
tags:
篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了SciPy笔记（有初始条件的常微分方程）相关的知识，希望对你有一定的参考价值。
参考技术A 3-1有初始条件的常微分方程import matplotlib.pyplot as plt
from IPython.display import Image
from scipy import *
import scipy.linalg as la
from scipy.integrate import odeint, ode
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
def f(t,y):
return y*np.sin(t)
def analytical_solution(t):
return np.e**(1-np.cos(t))
result=integrate.solve_ivp(f,(0,5),[1],dense_output=True)
#integrate.solve_ivp(f,(t_0,t_1),y_0,method='RK45',dense_output=False)
#(t_0, t_1)为参数t的范围，y0为初值条件。
#method='RK45'为选择的方法，默认为5（4）阶龙格库塔法。
#dense_output选择的是是否稠密输出,默认为否。
t=np.linspace(0,5,1001)
plt.plot(t,analytical_solution(t),label='解析解',color='g')
plt.plot(t,result.sol(t)[0],label='数值解',color='m')
plt.grid()
plt.legend()
plt.xlabel('t')
plt.ylabel('y(t)')
plt.title('有初始条件的常微分方程',fontsize =25)
plt.show()
计算方法带初值的常微分方程解法
实现了一些常见的带初值的常微分方程的算法。
/// <summary> /// 欧拉迭代公式,求出区间(x0,x]中每隔步长h的精度为e的近似值 /// </summary> /// <param name="fun">fun为x的函数即 dy/dx=fun(x,y)</param> /// <param name="N">将区间分为N段</param> /// <param name="x0">f(x0)=y0</param> /// <param name="y0">f(x0)=y0</param> /// <param name="x">所求区间的右端点</param> /// <param name="e">迭代精度</param> /// <returns>返回区间每隔h的函数值</returns> public static double[] Euler_e(Func<double, double, double> fun, int N, double x0, double y0, double x, double e = 1e-10) { double[] res = new double[N]; double h = (x - x0) / N; for (int i = 0; i < N; i++) { x0 = x0 + h; double yn = y0; if (i - 1 >= 0) yn = res[i - 1]; res[i] = yn + h * fun(x0 - h, yn); double res_e = 0; while (Math.Abs(res_e - res[i]) >= e) { res_e = res[i]; res[i] = yn + h / 2 * (fun(x0 - h, yn) + fun(x0, res_e)); } } return res; } /// <summary> /// 欧拉迭代公式,求出区间(x0,x]中每隔步长h的精度为e的近似值 /// </summary> /// <param name="fun">fun为x的函数即 dy/dx=fun(x)</param> /// <param name="N">将区间分为N段</param> /// <param name="x0">f(x0)=y0</param> /// <param name="y0">f(x0)=y0</param> /// <param name="x">所求区间的右端点</param> /// <param name="e">迭代精度</param> /// <returns>返回区间每隔h的函数值</returns> public static double[] Euler_e(Func<double, double> fun, int N, double x0, double y0, double x, double e = 1e-10) { double[] res = new double[N]; double h = (x - x0) / N; for (int i = 0; i < N; i++) { x0 = x0 + h; double yn = y0; if (i - 1 >= 0) yn = res[i - 1]; res[i] = yn + h * fun(x0 - h); double res_e = 0; while (Math.Abs(res_e - res[i]) >= e) { res_e = res[i]; res[i] = yn + h / 2 * (fun(x0 - h) + fun(x0)); } } return res; } //欧拉预估-矫正公式 //求出区间(x0,x]中每隔步长h的近似值 //fun为x的函数即 dy/dx=fun(x,y) //f(x0)=y0 public static double[] Euler_pre(Func<double, double, double> fun, int N, double x0, double y0, double x) { double[] res = new double[N]; double h = (x - x0) / N; for (int i = 0; i < N; i++) { x0 = x0 + h; double yn = y0; if (i - 1 >= 0) yn = res[i - 1]; res[i] = yn + h * fun(x0 - h, yn); res[i] = yn + h / 2 * (fun(x0 - h, yn) + fun(x0, res[i])); } return res; } //欧拉预估-矫正公式 //求出区间(x0,x]中每隔步长h的近似值 //fun为x的函数即 dy/dx=fun(x) //f(x0)=y0 public static double[] Euler_pre(Func<double, double> fun, int N, double x0, double y0, double x) { double[] res = new double[N]; double h = (x - x0) / N; for (int i = 0; i < N; i++) { x0 = x0 + h; double yn = y0; if (i - 1 >= 0) yn = res[i - 1]; res[i] = yn + h * fun(x0 - h); res[i] = yn + h / 2 * (fun(x0 - h) + fun(x0)); } return res; } //二阶龙格-库塔算法 //求出区间(x0,x]中每隔步长h的近似值 //fun为x的函数即 dy/dx=fun(x,y) //f(x0)=y0 public static double[] R_K2(Func<double, double, double> fun, int N, double x0, double y0, double x) { double[] res = new double[N]; double h = (x - x0) / N; for (int i = 0; i < N; i++) { double yn = y0; if (i - 1 >= 0) yn = res[i - 1]; double k1 = h * fun(x0, yn); double k2 = h * fun(x0 + 2.0 / 3 * h, yn + 2.0 / 3 * k1); res[i] = yn + 1.0 / 4 * (k1 + 3 * k2); x0 += h; } return res; } //二阶龙格-库塔算法 //求出区间(x0,x]中每隔步长h的近似值 //fun为x的函数即 dy/dx=fun(x) //f(x0)=y0 public static double[] R_K2(Func<double, double> fun, int N, double x0, double y0, double x) { double[] res = new double[N]; double h = (x - x0) / N; for (int i = 0; i < N; i++) { double yn = y0; if (i - 1 >= 0) yn = res[i - 1]; double k1 = h * fun(x0); double k2 = h * fun(x0 + 2.0 / 3 * h); res[i] = yn + 1.0 / 4 * (k1 + 3 * k2); x0 += h; } return res; } //四阶龙格-库塔算法 public static double[] R_K4(Func<double, double, double> fun, int N, double x0, double y0, double x) { double[] res = new double[N]; double h = (x - x0) / N; for (int i = 0; i < N; i++) { double yn = y0; if (i - 1 >= 0) yn = res[i - 1]; double k1 = h * fun(x0, yn); double k2 = h * fun(x0 + 0.5 * h, yn + 0.5 * k1); double k3 = h * fun(x0 + 0.5 * h, yn + 0.5 * k2); double k4 = h * fun(x0 + h, yn + k3); res[i] = yn + 1.0 / 6 * (k1 + 2 * (k2 + k3) + k4); x0 += h; } return res; } //四阶龙格-库塔算法re public static double[] R_K4(Func<double, double> fun, int N, double x0, double y0, double x) { double[] res = new double[N]; double h = (x - x0) / N; for (int i = 0; i < N; i++) { double yn = y0; if (i - 1 >= 0) yn = res[i - 1]; double k1 = h * fun(x0); double k2 = h * fun(x0 + 0.5 * h); double k3 = h * fun(x0 + 0.5 * h); double k4 = h * fun(x0 + h); res[i] = yn + 1.0 / 6 * (k1 + 2 * (k2 + k3) + k4); x0 += h; } return res; }
以上是关于SciPy笔记（有初始条件的常微分方程）的主要内容，如果未能解决你的问题，请参考以下文章