常微分方程组(边值)

合集下载

重要：MATLAB常微分方程(组)数值解法

Matlab常微分方程求解问题分类
边值问题：
初值问题：
• 定解附加条件在自变量的一端
• 一般形式为： y' f (x, y)
y(a)
y0
• 初值问题的数值解法一般采用步进法，如 Runge-Kutta法
➢ 在自变量两端均给定附加条件
y' f (x, y)
➢ 一般形式：y(a)y1, y(b)y2
1.根据常微分方程要求的求解精度与速度要求
求解初值问题：
y
'
y
2x y
y ( 0 ) 1
(0x1)
比较ode45和ode23的求解精度和速度
ode45和ode23的比较-1
function xODE clear all clc
format long
y0 = 1; [x1,y1] = ode45(@f,[0,1],y0); [x2,y2] = ode23(@f,[0,1],y0); plot(x1,y1,'k-',x2,y2,'b--') xlabel('x') ylabel('y')
rD = k(3)*C(2)-k(5)*C(4);
rE = k(4)*C(3)+k(5)*C(4);
% Mass balances dCdt = [rA; rB; rC; rD; rE];
三个串联的CSTR等温反应器（例4-3）
function IsothermCSTRs clear all clc CA0 = 1.8; % kmol/m^3 CA10 = 0.4; % kmol/m^3 CA20 = 0.2; % kmol/m^3 CA30 = 0.1; % kmol/m^3 k = 0.5; % 1/min tau = 2; stoptime = 2.9; % min [t,y] = ode45(@Equations,[0 stoptime],[CA10 CA20 CA30],[],k,CA0,tau); disp(' Results:') disp(' t CA1 CA2 CA3') disp([t,y]) plot(t,y(:,1),'k--',t,y(:,2),'b:',t,y(:,3),'r-') legend('CA_1','CA_2','CA_3') xlabel('Time (min)') ylabel('Concentration') % -----------------------------------------------------------------function dydt = Equations(t,y,k,CA0,tau) CA1 = y(1); CA2 = y(2); CA3 = y(3); dCA1dt = (CA0-CA1)/tau - k*CA1; dCA2dt = (CA1-CA2)/tau - k*CA2; dCA3dt = (CA2-CA3)/tau - k*CA3; dydt = [dCA1dt; dCA2dt; dCA3dt];

边值问题的数值解法

估计式
M b a 2 y xk y k h ，k 1， 2，，n 1。 96
2
y 4 x 。因此，当 h 0 时，差分方程的解收敛到微分方其中 M max a x b
y f x，y，y， y x，y sk，
这里的 s k 为
（8.6.3）
y
在处的斜率。令 z y ，上述二阶方程可降为一阶方程组
y z， z f x，y，z ，
（8.6.4）
y a ，z a sk。
计算结果表明打靶法的效果是很好的，计算精度取决于所选取的初值问题数
值方法的阶和所选取的步长 h 的大小。不过，打靶法过分依赖于经验，选取试射值，有一定的局限性。
第八章常微分方程数值解法
8.6.2 差分方法
差分方法是解边值问题的一种基本方法，它利用差商代替导数，将微分方程离散化为线性或非线性方程组（即差分方程）来求解。先考虑线性边值问题（8.6.2）的差分法。将区间 a，b 分成 n 等分，子区间的
s2
，同理得到 yb，s2 ，再判断它是否满足精度要求
y b，s2 。如此重复，直到某个 s 满足 y b，sk ，此时得到 k
的 y xi 和 yi z xi 就是边值问题的解函数值和它的一阶导数值。上述方程好比打靶， s k 作为斜率为子弹的发射，y b 为靶心，故称为打靶法。
y xy 4 y 12 x 2 3x， 0 x 1， y 0 0，y 1 2，
其解的解析表达式为 y
x x 4 x 。来自解先将该线性边值问题转化为两个初值问题
xy1 4 y1 12 x 2 3 x， y1 1 0， y1 0 0，y1 xy2 4 y2 0， y2 1 1。 y2 0 0，y2

解常微分方程组-边界值问题(Boundary

第十章解常微分方程組-邊界值問題(Boundary Value Problems for ODE)本章探討的邊界值問題模型如下:x''(t) = f(t, x(t), x'(t))x(a) = α , x(b) = β .考慮當函數f(t, x(t), x'(t))是線性(或非線性)時, 求解之方法:1. Discretization Method( finite difference approximations)2. Shooting Method在本章中包含Matlab 的m-file1. finitedf.m2.shoot.m3. shootnl.m將須要的m-file之檔案夾加入搜尋路徑中path('d:\numerical', path)註: 如果你有安裝Matlab Notebook 要執行下列input cells (綠色指令敘述)之前必須先執行上面的cell –[path (…) ]藍色的內容是Matlab [output cells]1. finitedf.m –利用finite difference approximationsx'(t) ~ (x(t+h) - x(t-h) ) / (2 h) ;x''(t) ~ (x(t+h) - 2 (t) + x(t-h) )/ h2 ;來估算邊界值問題x''(t) = f(t, x(t), x'(t))= u(t) t + v(t) x(t) + w(t) x'(t)x(a) = α , x(b) = β. 的數值解 .下列的m-file(finitedf.m)已經把例題中的u(t), v(t), w(t)定義了, 使用者套用此函數時, 記得更正為自己需要的.type finitedf.mfunction rs= finitedf(x0,xn, ta,tb,n)%to solve x''(t) = f(t, x(t), x'(t))=u(t) t + v(t) x(t) + w(t) x'(t)%by finite difference approximations, x0 and xn are the boundary values %t is in[ta,tb].ut = inline('exp(t) - 3*sin(t)', 't') ;vt = -1 ;wt = 1 ;h=(tb - ta)/n ;a=zeros(1,n); b=a; c=a; d=a;for i=1:n-1t = ta + i*h ;a(i) = -( 1+ h * wt/2) ;d(i) = 2 + h^2 * vt;c(i) = -( 1- h * wt/2) ;b(i) = -h^2* ut(t) ;end ;b(1) = b(1) - a(1)* x0 ;b(n-1) = b(n-1) - c(n-1)* xn ;for i=1:n-1a(i) = a(i+1) ;end ;y = trigauss(a, d, c,b,n-1) ;rs = y ;例題 1: 解邊界值問題x''(t) = e t - 3 sin(t) + x'(t) - x(t) ,x(1) = 1.09737491 , x(2) = 8.63749661ta = 1 ; tb = 2; n = 99;h = (tb- ta)/n ; error = zeros(1,n) ;x0 = 1.09737491 ; xn = 8.63749661 ;y = finitedf(x0, xn, ta, tb, n);for i = 1 : n-1t = ta + i * h ;error(i) = exp(t)- 3*cos(t)- y(i) ;end ;format longfprintf('\t t \t x(t) \t error\n') disp([1 x0 0])for i = 9 :9:n-1t = ta + i * h ;disp([t y(i) error(i)])end ;disp([2 xn 0])t x(t) error1.00000000000000 1.09737491000000 01.09090909090909 1.59194209918682 -0.000000365918851.181818181818182.12256804212776 -0.000000694314561.272727272727272.68955333563720 -0.000000963621881.36363636363636 3.29334396781902 -0.000001156328151.45454545454545 3.93457004989182 -0.000001259030391.54545454545455 4.61408706401417 -0.000001263391911.63636363636364 5.33301967061334 -0.000001167139561.72727272727273 6.09280813414278 -0.000000975096821.81818181818182 6.89525744460769 -0.000000700247811.90909090909091 7.74258923378422 -0.000000364826842.00000000000000 8.63749661000000 02. shoot.m –利用shooting method 來估算邊界值問題此方法是將邊界值問題轉為初值問題:x''(t) = f(t, x(t), x'(t)) = u(t) t + v(t) x(t) + w(t) x'(t)x(a) = α , x'(a) = z.其解x(t)在b的數值x(b)可視為z的函數φ(z) .希望好的z值能使得φ(z) =β.當φ(z)是線性時, 考慮函數g(t)=λx1(t) + (1-λ) x2(t) ,其中x1(t)與x2(t)分別是x'(a) = z1與.x'(a) = z2初值問題的解.解此初值問題可利用前面的函數rk4sys.mtype shoot.mfunction rs= shoot(x0,xn, ta,tb,n)%to solve x''(t) = f(t,x(t),x'(t))=u(t) t+ v(t) x(t)+ w(t) x'(t) %by finite difference approximations, x0 and xn are the%boundary values, t is in[ta,tb].ut = inline('exp(t) - 3*sin(t)', 't') ;vt = -1 ;wt = 1 ;h=(tb - ta)/n ; m=5 ;x=[1 x0 0 x0 1]';xt=zeros(m,n);xtnew=zeros(1,n) ;xt=rk4sysnew(x,ta,tb,n) ;xb1=xt(2,n) ;xb2=xt(4,n) ;p = (xn-xb2) / (xb1-xb2) ;q = 1-p ;xtnew=p*xt(2,1:n) +q*xt(4,1:n) ;%return the approximation solutionrs = xtnew ;type fxsys.mfunction fx=fxsys(t,X)%compute values of function F(t,X)%F(t,X) is defined by user%In this example, variable X includes t which is%the first component X(1)fx=zeros(length(X),1);fx(1) = 1 ;fx(2) = X(3);fx(3) = exp(X(1)) - 3*sin(X(1)) + X(3) -X(2);fx(4) = X(5) ;fx(5) = exp(X(1)) - 3*sin(X(1)) + X(5) -X(4); ;例題 2:利用shooting method解邊界值問題x''(t) = e t - 3 sin(t) + x'(t) - x(t) ,x(1) = 1.09737491 , x(2) = 8.63749661ta = 1 ; tb = 2; n = 99;h = (tb- ta)/n ; error = zeros(1,n) ;x0 = 1.09737491 ; xn = 8.63749661 ;y = shoot(x0, xn, ta, tb, n);for i = 1 : n-1t = ta + i * h ;error(i) = exp(t)- 3*cos(t)- y(i) ;end ;format longfprintf('\t t \t x(t) \t error\n') disp([1 x0 0])for i = 9 :9:n-1t = ta + i * h ;disp([t y(i) error(i)])end ; disp([2 xn 0])t x(t) error1.00000000000000 1.09737491000000 01.09090909090909 1.59194173254025 0.000000000727731.181818181818182.12256734722770 0.000000000585511.272727272727272.68955237158771 0.000000000427611.36363636363636 3.29334281123716 0.000000000253711.45454545454545 3.93456879079791 0.000000000063531.54545454545455 4.61408580076543 -0.000000000143171.63636363636364 5.33301850384033 -0.000000000366551.72727272727273 6.09280715965267 -0.000000000606711.81818181818182 6.89525674522358 -0.000000000863701.90909090909091 7.74258887009486 -0.000000001137472.00000000000000 8.63749661000000 03. shootnl.m –利用shooting method 來估算非線性的邊界值問題x''(t) = f(t, x(t), x'(t)) , x(a)= α , x(b)=β同樣地, 將邊界值問題轉為初值問題:x(a) = α , x'(a) = z.其解x(t)在b的數值x(b)可視為z的函數φ(z) .希望好的z值能使得φ(z) =β , 採用secant method .type shootnl.mfunction rs= shootnl(x0,xn, ta,tb,n)%to solve nonlinear x''(t) = f(t,x(t),x'(t)) boundary-Value problem%by Shooting method, x0 and xn are the boundary values,%t is in[ta, tb].%%x''(t)= -x'(t)^2 -x(t) + ln(t), t is in [1,2], x(1)=0, x(2) = ln2 .%exact solution x(t) = ln(t) .% x1(t)=t, x1'(t)=1, x2(t)= x(t), x2'(t)= x3(t)% x3'(t)= x2''(t)= -x3(t)^2 - x2(t) +ln(t)h=(tb - ta)/n ; m=3 ;%for some zz=[0.5 1] ; pz=[0 0] ; %initialepi=10^(-10) ;for j=1: 2x=[1 x0 z(j)]' ; %initialsxt=zeros(m,n) ;xt=rk4sysnew(x,ta,tb,n) ;pz(j)=xn -xt(2,n) ;if pz(j)< epirs = xt(2, 1:n) ;return %end ;end ;%compute the next z by secant methoditer= 15; tol = 10^-5 ; err=0.0 ;i =3 ;while(i <= iter & err < tol)slp= (z(i-1)-z(i-2))/(pz(i-1)-pz(i-2)) ;nextz = z(i-1)- slp*pz(i-1);err = abs(nextz - pz(i-1)) ;z=[z nextz] ;x=[1 x0 nextz]' ; %initialsxt=zeros(m,n) ;xt=rk4sysnew(x,ta,tb,n);pz(i)=xn -xt(2,n) ;if (err < tol | pz(i)<epi )rs = xt(2, 1:n ) ;returnend ;i= i+1 ;endif i>iterprintf(' Secant method fails in Shooting method \n' ) endrs = xt(2, 1:n); %returnstype fxsys.mfunction fx=fxsys(t,X)%compute values of function F(t,X)%F(t,X) is defined by user%In this example, variable X includes t which is%the first component X(1)fx=zeros(length(X),1);fx(1) = 1 ;fx(2) = X(3);fx(3) = -X(3).^2 - X(2) +log(X(1)) ;例題 3:利用shooting method 來估算非線性的邊界值問題x''(t) = -x'(t)2 - x(t) + ln(t) 1≤ t ≤ 2x(1) = 0 , x(2) = ln2 .ta = 1 ; tb = 2; n = 10;h = (tb- ta)/n ; error = zeros(1,n) ;x0 = 0.0 ; xn = log(2) ;y = shootnl(x0, xn, ta, tb, n);for i = 1 : n-1t = ta + i * h ;error(i) = log(t)- y(i) ;end ;format longfprintf('\t t \t x(t) \t error\n') fprintf('\t 1 \t 0 \t 0\n') for i = 1 :n-1t = ta + i * h ;disp([t y(i) error(i)])end ; disp([2 xn 0])t x(t) error1 0 01.10000000000000 0.09531000323697 0.000000176567361.20000000000000 0.18232131442583 0.000000242368121.30000000000000 0.26236400845782 0.000000256009671.40000000000000 0.33647199127437 0.000000245346841.50000000000000 0.40546488415551 0.000000223952651.60000000000000 0.47000343073393 0.000000198511811.70000000000000 0.53062807876778 0.000000172294391.80000000000000 0.58778651805356 0.000000146848561.90000000000000 0.64185376332331 0.000000122849092.00000000000000 0.69314718055995 0。

两点边值问题方程

两点边值问题方程两点边值问题是一种求解微分方程的方法，它涉及到两个边界条件。

假设我们有一个一阶常微分方程dy/dx = f(x,y)，我们需要找到满足两个边界条件y(a) = alpha 和y(b) = beta 的解。

两点边值问题的解法通常包括以下步骤：1. 定义一个初始猜测值y0(x)。

2. 使用数值方法（如欧拉法、龙格-库塔法等）求解微分方程，得到新的解y1(x)。

3. 检查新的解是否满足边界条件。

如果满足，则找到了解；否则，返回步骤2，使用新的解作为初始猜测值继续求解。

下面是一个使用Python实现两点边值问题的示例代码：```pythonimport numpy as npfrom scipy.integrate import odeint# 定义微分方程dy/dx = f(x,y)def f(x, y):return x * y - 1# 定义两个边界条件y(a) = alpha 和y(b) = betaa, b, alpha, beta = 0, 1, 1, 0# 定义初始猜测值y0(x)y0 = np.array([0.5, 0.5])# 使用数值方法求解微分方程def solve_two_point_boundary_value_problem(f, a, b, alpha, beta, y0, tol=1e-6, max_iter=100): for i in range(max_iter):y = odeint(f, y0, [a, b])if np.allclose(y[:1], alpha) and np.allclose(y[-1], beta):return y[1:-1]y0 = y[1:-1]raise ValueError("Solution not found within the specified tolerance and maximum iterations.")# 求解两点边值问题solution = solve_two_point_boundary_value_problem(f, a, b, alpha, beta, y0)print("Solution:", solution)```在这个示例中，我们使用`odeint`函数求解微分方程，并使用`np.allclose`检查新的解是否满足边界条件。

常微分方程的求解

18—1 常微分方程数值解法2§1 引言§2 Euler 方法§3 Runge －Kutta 方法§4 单步法的收敛性与稳定性§5 线性多步法§6 方程组与高阶方程的情况§7 边值问题的数值解法3§1 引言微分方程：关于一个未知函数的方程，方程中含有未知函数的（偏）导数，以及自变量等，其中关于未知函数导数的最高次数称为微分方程的阶数.例如：0)()(')()(''=++−x c y x b y x a x y4实际中，很多问题的数学模型都是微分方程. 常微分方程作为微分方程的基本类型之一，在理论研究与工程实际上应用很广泛. 很多问题的数学模型都可以归结为常微分方程. 很多偏微分方程问题，也可以化为常微分方程问题来近似求解.微分方程的应用情况5对于一个常微分方程：'(,) ,[,]dy y f x y x a b dx==∈为了使解存在，一般要对函数f 施加限制条件，例如要求f 对y 满足Lipschitz 条件：1212(,)(,)f x y f x y L y y −≤−6同时，一个有解的微分方程通常会有无穷多个解例如cos() sin(),dyx y x a a R dx=⇒=+∀∈为了使解唯一，需要加入一个限定条件. 通常会在端点出给出，如下面的初值问题：(,),[,]()dyf x y x a b dx y a y ⎧=∈⎪⎨⎪=⎩7常微分方程的解是一个函数，但是，只有极少数特殊的方程才能求解出来，绝大多数是不可解的.并且计算机没有办法对函数进行运算. 一般考虑其近似解法，一种是近似解析法，如逼近法、级数解法等，另一种是本章介绍的数值解法.8§2 Euler 方法92－1 Euler 公式对常微分方程初值问题：⎩⎨⎧==00')(),(y x y y x f y 数值求解的关键在于消除其中的导数项——称为离散化. 利用差商近似逼近微分是离散化的一个基本途径.10现在假设求解节点为),,1,0(m i ih a x i "=+=，其中ma b h −=为步长，这些节点相应的函数值为)(,),(1m x y x y ". 在点n x 处，已知))(,()('n n n x y x f x y =用n x 的向前差商nn n n x x x y x y −−++11)()(近似代替)('n x y ，如§1，则得到所谓的Euler 公式1(,)n n n n y y hf x y +=+——单步、显式格式11Euler 公式的局部截断误差：假设)(n n x y y =情况下，11)(++−n n y x y 称为局部截断误差.'''2311''23()()()()()2()(,()(()))2n n n n n n n n n y x y x y y x hy x h O h y x h y x f x y x h O h ++−=+++−−=+故有)(2)(''211n n n x y h y x y ≈−++. 122－2 后退的Euler 公式同样对常微分方程初值问题，在1+n x 点，已知))(,()(111'+++=n n n x y x f x y ，如果用向后差商hx y x y n n )()(1−+代替)(1'+n x y ，则得到后退的Euler 公式：111(,)n n n n y y hf x y +++=+——单步、隐式格式13相对于以上可以直接计算1+n y 的Euler 公式（显式），上式是隐式公式. 一般来讲，显式容易计算，而隐式具有更好的稳定性.求解上述公式，通常使用迭代法：对于给定的初值)0(1+n y，计算(1)()111(,)(0,1,)k k n n n n y y f x y k ++++=+="，如果)(1lim k n k y +∞→收敛，则其极限必满足上述后退Euler 公式.14局部截断误差：假设)(n n x y y =，则),()(111++++=n n n n y x hf x y y .由于)]()[,())(,(),(1111111+++++++−+=n n n y n n n n x y y x f x y x f y x f η且''''2111(,())()()()()n n n n n f x y x y x y x hy x O h +++==++15则有'2''31111(,)[()]()()()()n y n n n n n n y hf x y y x y x hy x h y x O h η++++=−++++将此式减去式2'''31()()()()()2n n n n h y x y x hy x y x O h +=+++ 可得，2''311111()(,)[()]()()2n n y n n n n h y x y hf x y x y y x O h η+++++−=−−+16考虑到21111(,)()1(,)y n y n hf x O h hf x ηη++=++−，则有22''3''11()()()()22n n n n h h y x y y x O h y x ++−=−+≈−172－3 梯形公式由于上述两个公式的局部截断误差绝对值相等，符号相反，故求其算术平均得到梯形公式：111[(,)(,)]2n n n n n n hy y f x y f x y +++=++——单步、隐式格式18梯形法同样是隐式公式，可用下列迭代公式求解：(0)1(1)()111(,)[(,)(,)]2n n n n k k n n n n n n y y hf x y h y y f x y f x y +++++⎧=+⎪⎨=++⎪⎩局部截断误差：类似于后退Euler ，可计算出)(12)('''311n n n x y h y x y −≈−++192－4 改进的Euler 公式上述用迭代法求解梯形公式虽然提高了精度，但计算量也很大. 实际上常采用的方法是，用Euler 公式求得初始值（预测），然后迭代法仅施行一次（校正）——改进的Euler 公式：1111(,)[(,)(,)]2n n n n n n n n n n y y f x y hy y f x y f x y ++++⎧=+⎪⎨=++⎪⎩20估计上式中第二式当1+n y 为准确值时的局部截断误差：''11113(3)()()(()[()()])2()12n n n n n n n hy x y y x y x y x y x hy x ++++−=−++≈−212－5 Euler 两步公式如果用中心差商hx y x y n n 2)()(11−+−代替)('n x y ，则得Euler 两步公式112(,)n n n n y y hf x y +−=+——两步、显式格式22假设1−n y 及n y 均为准确值，利用Taylor 展式容易计算Euler 两步公式的局部截断误差为：11113(3)()()(()2(,()))()3n n n n n n n y x y y x y x hf x y x h y x +++−−=−+≈23此式与梯形公式相结合，得到如下的预测－校正公式：111112(,)[(,)(,)]2n n n n n n n n n n y y hf x y hy y f x y f x y −++++⎧=+⎪⎨=++⎪⎩假设第一式中的1−n y 及n y ，以及第二式中的n y 及1+n y 均是准确值，则有，2441)()(1111−≈−−++++n n n n y x y y x y 从而可得以下的事后估计式，111111114()()51()()5n n n n n n n n y x y y y y x y y y ++++++++⎧−≈−−⎪⎪⎨⎪−≈−⎪⎩25可以期望，以上式估计的误差作为计算结果的补偿，可以提高计算精度.以n p 及n c 分别表示第n 步的预测值和校正值，则有以下的“预测－改进－校正－改进”方案（其中在1+n p 与1+n c 尚未计算出来的前提下，以n n c p −代替11++−n n c p ：26预测：'112n n n hy y p +=−+预测的改进：)(5411n n n n c p p m −−=++计算：),(11'1+++=n n n m x f m校正：)(2'1'1++++=n n n n m y hy c校正的改进：)(511111++++−+=n n n n c p c y计算：),(11'1+++=n n n y x f y27例用Euler 方法求解初值问题2'[0,0.6](0)1y y xy x y ⎧=−−∈⎨=⎩取0.2h =，要求保留六位小数. 解：Euler 迭代格式为2210.2()0.80.2k k k k k k k k y y y x y y x y +=+−−=−因此2821000(0.2)0.80.20.8y y y x y ≈=−= 22111(0.4)0.80.20.6144y y y x y ≈=−=23222(0.6)0.80.20.461321y y y x y ≈=−=29例用改进的Euler 方法求解初值问题2'sin 0[0,0.6](0)1y y y x x y ⎧++=∈⎨=⎩取0.2h =，求(0.2),(0.4)y y 的近似值，要求保留六位小数.解：改进的Euler 格式为212211110.2(sin )0.2(sin sin )2k k k k k k k k k k k k k y y y y x y y y y x y y x +++++⎧=+−−⎪⎨=+−−−−⎪⎩30即，222110.820.08sin 0.1(0.80.2sin )sin k k k k k k k k y y y x y y x x ++=−−−则有1(0.2)0.807285y y ≈=，2(0.4)0.636650y y ≈=31§3 Runge －Kutta 方法Def.1如果一种方法的局部截断误差为)(1+p h O ，则称该方法具有p 阶精度. 323－2 Runge —Kutta 方法的基本思想上述的Taylor 级数法虽然可得到较高精度的近似公式，但计算导数比较麻烦. 这里介绍不用计算导数的方法.))(,()()()('1h x y h x f h x y hx y x y n n n n n θθθ++=+=−+——平均斜率.33如果粗略地以),(n n y x f 作为平均斜率，则得Euler 公式；如果以221K K +作为平均斜率，其中),(1n n y x f K =，),(112hK y x f K n n +=+，则得改进的Euler 公式.343－3 二阶的Runge －Kutta 方法对点n x 和)10(≤<+=+p ph x x n p n ，用这两点斜率的线性组合近似代替平均斜率，则得计算公式：11122121()(,)(,)n n n n n p n y y h K K K f x y K f x y phK λλ++⎧=++⎪=⎨⎪=+⎩35现确定系数p ,,21λλ，使得公式具有二阶精度. 因为，取n y 为()n y x ，则'1(,)(,())'()n n n n n nK f x y f x y x y x y === 再把2K 在),(n n y x 处展开，有36'21(,)(,)n p n n n n K f x y phK f x ph y phy +=+=++代入可得，'2''31122()()n n n n y y hy ph y O h λλλ+=++++'2(,)(,)(,)()n n x n n y n n n f x y f x y ph f x y phy O h =+⋅+⋅+'2(')(,)()n x y n n y ph f f y x y O h =+⋅+⋅+'''2()n n y ph y O h =+⋅+37相比较二阶Taylor 展开''2'12n n n n y h hy y y ++=+，有，⎪⎩⎪⎨⎧==+211221p λλλ满足此条件的公式称为二阶Runge －Kutta 公式.38可以验证改进的Euler 公式属于二阶Runge －Kutta 公式. 下列变形的Euler 公式也是二阶Runge －Kutta 公式：12121(,)(,)22n n n n n n y y hK K f x y h h K f x y K +⎧⎪=+⎪=⎨⎪⎪=++⎩393－4 三阶Runge －Kutta 公式同二阶Runge －Kutta 公式，考虑三点,,(01)n n p n q x x x p q ++≤≤≤试图用它们的斜率321,,K K K 的线性组合近似代替平均斜率，即有如下形式的公式：1112233121312()(,)(,)(,())n n n n n n n n y y h K K K K f x y K f x ph y phK K f x qh y qh rK sK λλλ+=+++⎧⎪=⎪⎨=++⎪⎪=+++⎩40把32,K K 在),(n n y x 处展开，通过与)(1+n x y 在n x 的直接Taylor 展式比较，可确定系数s r q p ,,,,,,321λλλ，满足下式，从而使得上述公式具有三阶精度，41特别地，2,1,1,21,32,61231=−======s r q p λλλ是其一特例.123232223311213161p q p q pqs r s λλλλλλλλ++=⎧⎪⎪+=⎪⎪⎪+=⎨⎪⎪=⎪⎪+=⎪⎩423－5 四阶Runge －Kutta 公式相同的方法，可以导出下列经典的四阶Runge －Kutta 公式：112341213243(22)6(,)(,)22(,)22(,)n n n n n n n n n n h y y K K K K K f x y h h K f x y K h h K f x y K K f x h y hK +⎧=++++⎪⎪=⎪⎪⎪=++⎨⎪⎪=++⎪⎪=++⎪⎩43例用经典四阶Runge —Kutta 方法求解初值问题'83[0,0.4](0)1y y x y =−⎧∈⎨=⎩，取0.2h =，求(0.4)y 的近似值，要求保留六位小数.解：四阶Runge —Kutta 格式为44112341211123122241330.2(22)6(,)830.2(,)83(0.1) 5.6 2.120.2(,)83(0.1) 6.32 2.372(,0.2)83(0.2) 4.208 1.578k k k k k k k k k k k kk k k k ky y K K K K K f x y y K f x y K y K yK f x y K y K y K f x y K y K y ++++⎧=++++⎪⎪==−⎪⎪⎪=+=−+=−⎨⎪⎪=+=−+=−⎪⎪⎪=+=−+=−⎩则10.5494 1.2016k k y y +=+，45故12(0.2) 2.3004,(0.4) 2.4654y y y y ≈=≈=.注：由准确解382()33xy x e −=−可得(0.2) 2.300792,(0.4) 2.465871y y ==46§5 线性多步法基本思想：在计算1+i y 之前，已计算出一系列的近似值i y y ,,1"，如果充分利用这些已知信息，可以期望会获得更高精度的)(1+i x y 的近似值1+i y .基本方法：基于数值积分与基于Taylor 展开的构造方法.475－1 基于数值积分的构造方法对方程),('y x f y =两边从i x 到1+i x 积分，则得∫++=+1),()()(1i ix x i i dxy x f x y x y 设)(x P r 是f (x , y )的插值多项式，由此可得以下的一般形式的计算公式：∫++=+1)(1i ix x r i i dxx P y y 48例取线性插值))(,())(,()(11111+++++−−+−−=i i i i ii i i i i r x y x f x x x x x y x f x x x x x P ，则得到梯形法：)],(),([2111+++++=i i i i i i y x f y x f hy y495－2 Adams 显式公式在区间],[1+i i x x 上利用r ＋1个数据点),(,),,(),,(11r i r i i i i i f x f x f x −−−−"构造插值多项式)(x P r ，由牛顿后插公式（注意到：j i j i j f f −Δ=∇）j i jrj j i r f j t th x P −=Δ⎟⎟⎠⎞⎜⎜⎝⎛−−=+∑0)1()(其中!)1()1(j j s s s j s +−−=⎟⎟⎠⎞⎜⎜⎝⎛". 50可得10rj i i rj i jj y y h f αΔ+−==+∑——Adams 显式公式其中1(1)j j t dt j α−⎛⎞=−⎜⎟⎝⎠∫，它可写成：∑=−++=rj ji rj i i f h y y 01β515－3 Adams 隐式公式在区间],[1+i i x x 上利用r ＋1个数据点),(,),,(),,(1111+−+−++r i r i i i i i f x f x f x "构造插值多项式)(x P r ，由牛顿后插公式101)1()(+−=+Δ⎟⎟⎠⎞⎜⎜⎝⎛−−=+∑j i jrj ji r f j t th x P 可得*11rj i i rj i j j y y h f α+−+==+Δ∑——Adams 隐式公式52其中01(1)jj t dt j −−⎛⎞α=−⎜⎟⎝⎠∫，它又可写成： *11ri i rj i j j y y h f β+−+==+∑535－4 Adams 预测－校正公式以r ＝3时的Adams 显式与隐式公式为例. 此时，显式公式为)9375955(243211−−−+−+−+=i i i i i i f f f f hy y 利用Taylor 展式，容易计算局部截断误差为)(720251)5(5i x y h . 54)5199(242111−−+++−++=i i i i i i f f f f hy y 同样利用Taylor 展开可得，其局部截断误差为5(5)19()720i h y x −. 隐式公式为55⎪⎩⎪⎨⎧+−++=−+−+=−−+++−−−+)519),(9(24)9375955(24211113211i i i i i i i i i i i i i f f f y x f hy y f f f f h y y 注利用2－5节的相同作法同样可以构造更精确的计算过程.可构造利用显式预测，隐式校正的计算公式：56§6 方程组与高阶方程的情形6－1 一阶方程组常微分方程初值问题为⎩⎨⎧==00)(),('y x y y x f y 此时T m y y y ),,(1"=，Tm f f f ),,(1"=. 此时上述的一切方法均可使用，只是注意y 与f 此时为向量.576－2 化高阶方程为一阶方程组解下列的m 阶方程()(1)'(1)(1)000000(,,',,)(),'(),,()m m m m y f x y y y y x y y x y yx y −−−⎧=⎨===⎩""令)1(21,,',−===m m y y y y y y "，则有58'12'23'1'12(,,,,)m m m m y y y y y yy f x y y y −⎧=⎪=⎪⎪⎨⎪=⎪⎪=⎩#"初始条件为：)1(00'002001)(,,)(,)(−===m m y x y y x y y x y "。

课件：级第四章 2 边值问题

y(a)
(a x b)
y(b)
例 3：传热问题建立微分方程：
d 2T 1 dT f (r) dr2 r dr
对流传热
建立边界条件：
a1T(b) b1T(b) T1
r=b
T1
●第三类边界条件
－给定边界处函数和导数共同满足的条件
●第三类边值问题
y f (x, y, y) (a x b)
●● ● ●●
y(x)
x =a
x =b
打靶法的几何说明
对于初值问题
y f x, y, y a x b
y(a)
y(a) m
m
m0
m1
mn
y(b) y(b)m0 y(b)m1 y(b)mn
y(b)m F(m)
合适的 m 值应满足：
y(b)m
即： F(m)
化标准形式：f (m) F(m) 0
1T
2
解：第一步：明确需要确定哪些函数值 u0,u1,u2,,uN,uN1
将
Ti
ui1
2ui h2
ui1
代入离散化方程
h2 ui1 2ui ui1 k g(Zi )
u0 2u1 u2
u1 2u2 u3 uN 1 2uN
h2
k
h2 k
u N 1
g (Z1 )
g(Z2 ) h2
●第一类边值问题
y f (x, y, y) (a x b)
y(a) y(b)
例 2：传热问题
绝热 r=b
建立微分方程：
d 2T 1 dT f (r) dr2 r dr
建立边界条件：
T (b) 0
●第二类边界条件－给定边界处导函数满足的条件

边值问题的数值解法在具体求解常微分方程时-2022年学习资料

中南大学数学科学与计算技术学院-第八章常微分方程数值解法-=323-z2=-x32+4y2?-y20=0, 20=0。-取h=0.02,用经典R-K法分别求这两个方程组解yx和y2x的计算值y1:和-y2i,然后按 8.6.6得精确解-6=,t2=0.x-y21-的打靶法计算值》：，部分点上的计算值、精确值和误差列于表8 12。-版核防行：小人学影学烧

中南大学数学科学与计算技术学院-第八章常微分方程数值解法-值得指出的是，对于线性边值问题86.2，一个简单实用的方法是用解-析的思想，将它转化为两个初值问题：-y"+pxyi+qxy =fx-ya=a,ya=0: 「片+px5+gxy2=0,-ly2a=a,y2a=l。-求得这两个初值问题的解yx和y2x,若y2b≠0 容易验证-a高-8.6.6-为线性两点边值问题8.6.2的解。-例8.7用打靶法求解线性边值问题-版核防行小人学影学烧
中南大学数学科学与计算技术学院-第八章常微分方程数值解法-y”+y-4y=12x2-3x,0<x<1,-1 0=0,y1=2,-其解的解析表达式为yX=x4+x。-解先将该线性边值问题转化为两个初值问题-y0=0, 1=0,-y2+y%-4y2=0,-y20=0,y1=1。-令乙1=2=y?，将上述两个边值问题分别降为一方程组初值问题-31=-x31+4y1+12x2-3x,-y,0=0,z10=0,-版权防行：小人学影学烧
中南大学数学科学与计算技术学院-第八章常微分方程数值解法-表8-12-Xi-yu-y2i-yx-y-yl-0.2--0.002407991-0.204007989-0.2016000053-,0.2016000 00-0.53×10-8-0.4--0.006655031-0.432255024-0.425600008 -0.4256000000-0.80x108-0.6-0.019672413-0.709927571-0. 2960000830.7296000000-0.83×108-0.145529585-1.06407038 -1.2096000058-1.2096000000-0.58x108-0.475570149-1.524 28455-2.00000000002.0000000000-例8.8用打靶法求解线性边值问题-4y"+y =2x3+16,-y2=8,y3=35/3。-要求误差不超过0.5×106，其解析解是yx=x2+8/x。解对应于8.6.4的初值问题为-版凤防行：小人学数：学烧

打靶法

用某种离散化数值步骤求出常微分方程边值问题在离散点上的近似解的方法。

各种实际问题导出不同类型的边值问题。

较简单的有二阶常微分方程两点边值问题：求函数y=y(x)，x∈【α,b】，使它满足微分方程和边值条件式中ƒ、g1、g2为已知函数；α与b为两个给定的端点。

较一般地有一阶常微分方程组两点边值问题:求N个函数使其满足微分方程组和边值条件式中诸ƒn、g i是已知函数;r为给定的自然数。

有些问题因求解区间是无穷区间而被称作奇异边值问题，相应的边界条件变为对解在无穷远处渐近行为的限制,例如,要求y(x)在区间【0，)上平方可积或要求当x趋于无穷时，y(x)趋于某极限值。

还有些实际问题因要求解满足多个点上的条件而被称作多点边值问题。

近年来，对反映边界层现象的奇异摄动边值问题提出了一些新的数值解法。

此外，关于存在多个解的分歧现象数值解问题也引起人们的注意。

打靶法主要思路是:适当选择和调整初值条件,(选什么)求解一系列初值问题，使之逼近给定的边界条件。

如果将描述的曲线视作弹道，那么求解过程即不断调整试射条件使之达到预定的靶子,所以称作打靶法或试射法,此类方法的关键是设计选取初值的步骤。

对非线性边值问题可通过下列步骤求数值解：①计算初值问题的数值解y1。

若g(y1(b),y姈(b))=B,近似地满足，则y1即为所求；否则进行②。

②计算初值问题的数值解y2,若g(y2(b),y娦(b))=B近似地满足，则y2即为所求；否则令m＝3进行③。

③将g(y(b),y┡(b))视为y(α)的函数,用线性逆插值法调整初值，即计算然后进行④。

④计算初值问题的数值解y m并进行判定：若b点边值条件近似地满足,则y m即为所求;否则令m+1崊m转向③继续计算直到满意为止。

特别地，若微分方程是线性的，则打靶法变成线性组合法，即根据常微分方程理论适当选取初值可得到一组线性独立解，利用它们的线性组合导出边值问题的解。

例如线性方程边值问题的数值解可通过两个初值问题数值解来实现。

常微分方程的边值问题

常微分方程的边值问题
常微分方程的边值问题（也称为常微分方程的定边值问题）是求解一个微分方程在一个给定的时间段上的特定解的问题，其中方程的解需要满足一些给定的边界条件。

这些边界条件通常指定了方程在时间间隔的起点和终点处的值，或者其他一些特定的时刻或位置上的值。

例如，一个常见的常微分方程的边值问题是求解一个二阶常微分方程：
y''(t) = f(t, y(t))
其中，y(t) 是未知函数，f(t, y) 是一个已知的函数。

这个问题需要在给定的时间段 [a, b] 上求解，并且需要满足以下的边界条件：
y(a) = y_a
y(b) = y_b
这里，y_a 和 y_b 是给定的数值。

这些边界条件指定了方程在时间间隔的起点和终点处的值。

常微分方程的边值问题在物理学、工程学、经济学等领域中都有广泛的应用。

解决常微分方程的边值问题需要使用数值解法或者解析解法，其中数值解法通常更为实用，因为它可以通过计算机程序来求解。

二阶常微分方程边值问题数值方法

其中 p( x),q( x)为,r已( x知) 函数，则由常微分方程的理论知，通过
变量替换总可以消去方程中的项，不妨y设变换后的方程为
y( x) q( x) y( x) r( x)
y(a) ,
y(b)
则近似差分方程成离散差分方程为
yi 1
2 yi h2
yi 1
qi
yi
ri
其中 qi q( xi ), ri r( xi ), i 1,2, , n. y0 ,
第一边界问题：
y0 , yn1
(8.9)
第二边界问题：
y1 y0 h , yn1 yn h
(8.10)
第三边界问题：
y1 (1 0h) y0 1h,
(1 0h) yn1 yn 1h
(8.11)
若 f ( x, y,是y) 的y线, y性函数时，f 可写成
f (x, y, y) p(x) y( x) q( x) y(x) r( x)
以
y
为待定参数。
0
对第三类边界问题，仍可转化为考虑初值问题(8.5)，取
y0 ,
y0 1 0 y0 ，以 y为0 待定参数。
8.2 有限差分法
将区间[a,b]进行等分：
h
ba, n1
xi
a ih, i 0,1,
,n 1,
设在
x xi , i 0,1, , n 1处的数值解为。 yi 用中心差分近似微分，即
而且还有误差估
计：
Ri
y( xi )
yi
M 24
h2
(
xi
a)(b xi )
其中 M max y(4。) ( x)
x[a ,b]

Matlab求解常微分方程边值问题的方法

Matlab 求解常微分方程边值问题的方法：bvp4c 函数常微分方程的边值问题，即boundary value problems ，简称BVP 问题，是指表达形式为(,)((),())0'=⎧⎨=⎩y f x y g y a y b 或(,,)((),(),)0'=⎧⎨=⎩y f x y p g y a y b p 的方程组（p 是未知参数），在MATLAB 中使用积分器bvp4c 来求解。

[命令函数]bvp4c[调用格式]sol=bvp4c(odefun,bcfun,solinit,options,p1,p2,…)sol 为一结构体，sol.x 、sol.y 、sol.yp 分别是所选择的网格点及其对应的y(x)与y'(x)数值; bvp4c 为带边值条件常微分方程积分器的函数命令；odefun 为描述微分方程组的函数文件；bcfun 为计算边界条件g(f(a),f(b),p)=0的函数文件；solinit 为一结构体，solinit.x 与solinit.y 分别是初始网格的有序节点与初始估计值，边界值条件分别对应a=solinit.x(l)和b=solinit.x(end)； options 为bvpset 命令设定的可选函数，可采用系统默认值；p1, p2…为未知参数。

例求常微分方程0''+=y y 在(0)2=y 与(4)2=-y 时的数值解。

[解题过程] 仍使用常用方法改变方程的形式：令1=y y ，21'=y y ，则原方程等价于标准形式的方程组1221⎧'=⎪⎨'=-⎪⎩y y y y ；将其写为函数文件twoode.m ；同时写出边界条件函数对应文件twobc.m ；分别使用结构solinit 和命令bvp4c 确定y-x 的关系；作出y-x 的关系曲线图。

[算例代码]solinit =bvpinit(linspace(0,4,5),[1 0]); % linspace(0,4,5)为初始网格，[1,0]为初始估计值 sol=bvp4c(@twoode,@twobc,solinit);% twoode 与twobc 分别为微分方程与边界条件的函数，solinit 为结构x=linspace(0,4); %确定x 范围y=deval(sol,x); %确定y 范围plot(x,y(1,:)); %画出y-x 的图形%定义twoode 函数（下述代码另存为工作目录下的twoode.m 文件）function dydx= twoode(x,y) %微分方程函数的定义dydx =[y(2) -abs(y(1))];%定义twobc 函数（下述代码另存为工作目录下的twobc.m 文件）function res= twobc(ya,yb); %边界条件函数的定义res=[ya(1);yb(1)+2];[运行结果]。

常微分方程边值问题解法

常微分方程边值问题解法
常微分方程边值问题解法：
常微分方程边值问题是指在一定区间内，给定一个微分方程的初始条件和边界条件，求解微分方程的解在这个区间内满足这些条件的问题。

常见的边值问题有两种类型：Dirichlet边界条件和Neumann边界条件。

解决常微分方程边值问题的方法有很多种，下面介绍其中两种常用的方法：
1. 有限差分法：
有限差分法是利用差分近似替代微分，将微分方程转化为一组代数方程。

首先将区间离散化，将连续的函数转化为离散的数值，然后利用中心差分、前向差分或后向差分的方法，将微分方程变为代数方程组，最后利用线性代数的方法求解这个方程组。

2. 有限元法：
有限元法是将区间划分为若干个小的子区间，将微分方程转化为一组局部的代数方程组，然后将这些方程组组合成整个问题的全局方程组。

有限元法可以适用于更加复杂的边值问题，但是需要更多的计算量和更高的数学水平。

总之，常微分方程边值问题的解法有很多种，需要根据具体情况选择不同的方法。

第8章--常微分方程边值问题的数值解法

第8章常微分方程边值问题的数值解法8.1 引言推论若线性边值问题()()()()()(),,(),()y x p x y x q x y x f x a x b y a y b αβ'''=++≤≤⎧⎨==⎩ (8.1.2) 满足(1) (),()p x q x 和()f x 在[,]a b 上连续； (2) 在[,]a b 上, ()0q x >, 则边值问题(8.1.1)有唯一解。

求边值问题的近似解，有三类基本方法：(1) 差分法(difference method)，也就是用差商代替微分方程及边界条件中的导数，最终化为代数方程求解；(2) 有限元法(finite element method)；(3) 把边值问题转化为初值问题，然后用求初值问题的方法求解。

8.2 差分法8.2.1 一类特殊类型二阶线性常微分方程的边值问题的差分法设二阶线性常微分方程的边值问题为(8.2.1)(8.2.2)()()()(),,(),(),y x q x y x f x a x b y a y b αβ''-=<<⎧⎨==⎩其中(),()q x f x 在[,]a b 上连续，且()0q x ≥.用差分法解微分方程边值问题的过程是：(i) 把求解区间[,]a b 分成若干个等距或不等距的小区间，称之为单元；(ii) 构造逼近微分方程边值问题的差分格式. 构造差分格式的方法有差分法, 积分插值法及变分插值法；本节采用差分法构造差分格式；(iii) 讨论差分解存在的唯一性、收敛性及稳定性；最后求解差分方程. 现在来建立相应于二阶线性常微分方程的边值问题(8.2.1), (8.2.2)的差分方程. ( i ) 把区间[,]I a b =N 等分，即得到区间[,]I a b =的一个网格剖分：011N N a x x x x b -=<<<<=,其中分点(0,1,,)i x a ih i N =+=，并称之为网格节点(grid nodes)；步长b a Nh -=. ( ii ) 将二阶常微分方程(8.2.2)在节点i x 处离散化：在内部节点(1,2,,1)i x i N =-处用数值微分公式2(4)1112()2()()()(),12i i i i i i i i y x y x y x h y x y x x h ξξ+---+''=-<< (8.2.3)代替方程(8.2.2)中()i y x ''，得112()2()()()()()()i i i i i i i y x y x y x q x y x f x R x h +--+-=+,(8.2.4) 其中2(4)()()12i i h R x y ξ=. 当h 充分小时，略去式(8.2.4)中的()i R x ，便得到方程(8.2.1)的近似方程1122i i i i i i y y y q y f h +--+-=,(8.2.5)其中(),()i i i i q q x f f x ==, 11,,i i i y y y +-分别是11(),(),()i i i y x y x y x +-的近似值, 称式(8.2.5)为差分方程(difference equation)，而()i R x 称为差分方程(8.2.5)逼近方程(8.2.2)的截断误差(truncation error). 边界条件(8.7.2)写成0,.N y y αβ==(8.2.6)于是方程(8.2.5), (8.2.6)合在一起就是关于1N +个未知量01,,,N y y y ，以及1N +个方程式的线性方程组：2211212211222111(2),(2),1,2,,1,(2).i i i i i N N N N q h y y h f y q h y y h f i N y q h y h f αβ-+----⎧-++=-⎪-++==-⎨⎪-+=-⎩(8.2.7)这个方程组就称为逼近边值问题(8.2.1), (8.2.2)的差分方程组(system of difference equations)或差分格式(difference scheme)，写成矩阵形式2211122222223332222222111(2)11(2)11(2)11(2)11(2)N N N N N N y q h h f y q h h f y q h h f y q h h f y q h h f αβ------⎡⎤⎡⎤-+-⎡⎤⎢⎥⎢⎥⎢⎥-+⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥-+=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥-+⎢⎥⎢⎥⎢⎥-+-⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦. (8.2.8)用第2章介绍的解三对角方程组的追赶法求解差分方程组(8.2.7)或(8.2.8), 其解01,,,N y y y 称为边值问题(8.2.1), (8.2.2)的差分解(difference solution). 由于(8.2.5)是用二阶中心差商代替方程(8.2.1)中的二阶微商得到的，所以也称式(8.2.7)为中心差分格式(centered-difference scheme).( iii ) 讨论差分方程组(8.2.7)或(8.2.8)的解是否收敛到边值问题(8.2.1), (8.2.2)的解，估计误差.对于差分方程组(8.2.7)，我们自然关心它是否有唯一解；此外，当网格无限加密，或当0h →时，差分解i y 是否收敛到微分方程的解()i y x . 为此介绍下列极值原理：定理8.2.1 (极值原理) 设01,,,N y y y 是给定的一组不全相等的数，设1122(),0,1,2,,i i i i i i i y y y l y q y q i N h +--+=-≥=.(8.2.9)(1) 若()0,1,2,,i l y i N ≥=, 则{}0Ni i y =中非负的最大值只能是0y 或N y ； (2) 若()0,1,2,,i l y i N ≤=, 则{}0Ni i y =中非正的最小值只能是0y 或N y .证只证(1) ()0i l y ≥的情形，而(2) ()0i l y ≤的情形可类似证明. 用反证法. 记{}0max i i NM y ≤≤=，假设0M ≥, 且在121,,,N y y y -中达到. 因为i y 不全相等，所以总可以找到某个00(11)i i N ≤≤-，使0i y M =，而01i y -和01i y +中至少有一个是小于M 的. 此时0000000011222()2.i i i i i i i i y y y l y q y h M M M q M q M h +--+=--+<-=-因为00,0i q M ≥≥，所以0()0i l y <, 这与假设矛盾，故M 只能是0y 或N y . 证毕！推论差分方程组(8.2.7)或(8.2.8)的解存在且唯一. 证明只要证明齐次方程组11202()0,0,1,2,,1,0,0i i i i i i i N y y y l y q y q i N h y y +--+⎧=-=≥=-⎪⎨⎪==⎩ (8.2.10)只有零解就可以了. 由定理8.7.1知，上述齐次方程组的解01,,,N y y y 的非负的最大值和非正的最小值只能是0y 或N y . 而00,0N y y ==，于是0,1,2,,.i y i N == 证毕！利用定理8.2.1还可以证明差分解的收敛性及误差估计. 这里只给出结果：定理8.2.2 设i y 是差分方程组(8.2.7)的解，而()i y x 是边值问题(8.2.1), (8.2.2)的解()y x 在i x 上的值，其中0,1,,i N =. 则有224()(),96i i i M h y x y b a ε=-≤-(8.2.11)其中(4)4max ()a x bM yx ≤≤=.显然当0h →时，()0i i i y x y ε=-→. 这表明当0h →时，差分方程组(8.2.7)或(8.2.8)的解收敛到原边值问题(8.7.1), (8.7.2)的解.例8.2.1 取步长0.1h =，用差分法解边值问题43,01,(0)(1)0,y y x x y y ''-=≤≤⎧⎨==⎩并将结果与精确解()()2222()3434x xy x e e ee x --=---进行比较.解因为110N h ==，()4,()3q x f x x ==, 由式(8.2.7)得差分格式221222112289(240.1)30.10.1,(240.1)30.1,2,3,,8,(240.1)30.10.9,i i i i y y y y y x i y y -+⎧-+⨯+=⨯⨯⎪-+⨯+=⨯=⎨⎪-+⨯=⨯⨯⎩0100y y ==, 00.1,1,2,,9i x ih i i =+==, 其结果列于表8.2.1.从表8.2.1可以看出, 差分方法的计算结果的精度还是比较高的. 若要得到更精确的数值解，可用缩小步长h 的方法来实现.8.2.2 一般二阶线性常微分方程边值问题的差分法对一般的二阶微分方程边值问题1212()()()()()(),,()(),()(),y x p x y x q x y x f x a x b y a y a y b y b αααβββ'''++=<<⎧⎪'+=⎨⎪'+=⎩ (8.2.12) 假定其解存在唯一.为求解的近似值，类似于前面的做法，( i ) 把区间[,]I a b =N 等分，即得到区间[,]I a b =的一个网格剖分：011N N a x x x x b -=<<<<=,其中分点(0,1,,)i x a ih i N =+=，步长b a Nh -=. ( ii ) 对式(8.2.12)中的二阶导数仍用数值微分公式2(4)1112()2()()()(),12i i i i i i i iy x y x y x h y x y x x h ξξ+---+''=-<<代替，而对一阶导数，为了保证略去的逼近误差为2()O h ，则用3点数值微分公式；另外为了保证内插，在2个端点所用的3点数值微分公式与内网格点所用的公式不同，即21112012000022212()()()(),,1,2,,1,263()4()()()(),,23()4()3()()(),.23i i i i i i i N N N N N N N N y x y x h y x y x x i N h y x y x y x h y x y x x h y x y x y x h y x y x x h ξξξξξξ+-----⎧-''''=-<<=-⎪⎪-+-⎪''''=+<<⎨⎪⎪-+''''=+<<⎪⎩(8.2.13) 略去误差，并用()i y x 的近似值i y 代替()i y x ，(),(),()i i i i i i p p x q q x f f x ===，便得到差分方程组1111221001221211(2)(),1,2,,1,2(34),2(43),2i i i i i i i i i N N N N p y y y y y q y f i N h hy y y y h y y y y h αααβββ-++---⎧-++-+==-⎪⎪⎪+-+-=⎨⎪⎪+-+=⎪⎩(8.2.14)其中(),(),(),1,2,,1i i i i i i q q x p p x f f x i N ====-, i y 是()i y x 的近似值. 整理得12021222211222121(23)42,(2)2(2)(2)2,1,2,,1,4(32)2.i i i i i i i N N N h y y y h hp y h q y hp y h f i N y y h y h αααααβββββ-+---+-=⎧⎪---++==-⎨⎪-++=⎩ (8.2.15)解差分方程组(8.2.15)，便得边值问题(8.2.12)的差分解01,,,N y y y .特别地, 若12121,0,1,0ααββ====，则式(8.2.12)中的边界条件是第一类边值条件：(),();y a y b αβ==此时方程组(7.7.16)为221112112211221211112(2)(2)2(2),(2)2(2)(2)2,2,3,,2,(2)2(2)2(2).i i i i i i i N N N N N N h q y hp y h f hp hp y h q y hp y h f i N hp y h q y h f hp αβ-+------⎧--++=--⎪---++==-⎨⎪---=-+⎩(8.2.16) 方程组(8.2.16)是三对角方程组，用第2章介绍的解三对角方程组的追赶法求解差分方程组(8.2.16)，便得边值问题(8.2.12)的差分解01,,,N y y y .( iii ) 讨论差分方程组(8.2.16)的解是否收敛到微分方程的解，估计误差. 这里就不再详细介绍.例8.2.2取步长/16h π=，用差分法求下列边值问题的近似解，并将结果与精确解进行比较.精确解是1()(sin 3cos )10y x x x =-+. 解因为(20)8N h π=-=，()1,()2,()cos p x q x f x x =-=-=, 由式(8.2.17)得差分格式()()()()()()()()()()()()()2122211222122216(2)216(1)216cos 16216(1)(0.3),216(1)2216(2)216(1)216cos 16,2,3,,6,216(1)2216(2)216cos 7i i i N N y yy y y i i y y πππππππππππππ-+--⎡⎤--⨯-++⨯-⎡⎤⎣⎦⎣⎦=--⨯-⨯-⎡⎤⎣⎦⎡⎤-⨯---⨯-++⨯-⎡⎤⎡⎤⎣⎦⎣⎦⎣⎦==⎡⎤-⨯---⨯-⎡⎤⎣⎦⎣⎦=()()16216(1)(0.1),ππ⎧⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪-+⨯-⨯-⎡⎤⎣⎦⎩080.3,0.1y y =-=-, 00.1,1,2,,9i x ih i i =+==, 其结果列于表8.2.2.8.3 有限元法有限元法(finite element method)是求解微分方程定解问题的有效方法之一，它特别适用在几何、物理上比较复杂的问题. 有限元法首先成功地应用于结构力学和固体力学，以后又应用于流体力学、物理学和其他工程科学. 为简明起见，本节以线性两点边值问题为例介绍有限元法.考虑线性两点边值问题()(8.3.1)(8.3.2)()()()()(),,(),(),Ly p x y x q x y x f x a x b y a y b αβ⎧''⎪=-+=≤≤⎨==⎪⎩其中1()0,()0,C [,]p x q x p a b >≥∈, ,C[,]q f a b ∈.此微分方程描述了长度为b a -的可变交叉截面(表示为()q x )的横梁在应力()p x 和()f x 下的偏差()y x .8.3.1 等价性定理记{}221C [,]()C [,],(),()a b y y y x a b y a y b αβ==∈==, 引进积分()22()()[()]()()2()()d baI y p x y x q x y x f x y x x '=+-⎰. (8.3.3)任取21()C [,]y y x a b =∈，就有一个积分值()I y 与之对应，因此()I y 是一个泛函(functional)，即函数的函数. 因为这里是,y y '的二次函数，因此称()I y 为二次泛函.对泛函(8.3.3)有如下变分问题(variation problem)：求函数21C [,]y a b ∈，使得对任意21C [,]y a b ∈, 均有()()I y I y ≥, (8.3.4)即()I y 在y 处达到极小, 并称y 为变分问题(8.3.4)的解.可以证明：定理8.3.1(等价性定理) y 是边值问题(8.3.1), (8.3.2)的解的充分必要条件是y 使泛函()I y 在21C [,]a b 上达到极小，即y 是变分问题(8.3.4)在21C [,]a b 上的解. 证 (充分性) 设21C [,]y a b ∈是变分问题()I y 的解；即y 使泛函()I y 在21C [,]a b 上达到极小，证明y 必是边值问题(8.3.1), (8.3.2)的解.设()x η是2C [,]a b 任意一个满足()()0a b ηη==的函数，则函数21()()()C [,]y x y x x a b αη=+∈，其中α为参数. 因为y 使得()I y 达到极小，所以()()I y I y αη+≥，即积分()22()()[()()]()[()()]2()[()()]baI y p x y x x q x y x x f x y x x dxαηαηαηαη''+=+++-+⎰作为α的函数，在0α=处取极小值()I y ，故d()0d I y ααηα=+=. (8.3.5) 计算上式，得()()()()()022(8.d()d d ()[()()]()[()()]2()[()()]d d 2()[()()]()2()[()()]()2()()d 2()()()()()()()()d .bab abaI y p x y x x q x y x x f x y x x x p x y x x x q x y x x x f x x x p x y x x q x y x x f x x x ααααηααηαηαηααηηαηηηηηη===+''=+++-+'''=+++-''=+-⎰⎰⎰3.6)利用分部积分法计算积分[][]()()()d ()()d ()()()()()()()d ()()()d ,bbaab ba abap x y x x x p x y x x p x y x x x p x y x x x p x y x x ηηηηη'''='''=-''=-⎰⎰⎰⎰代入式(8.3.6)，得()0(8.3.7)d()2()()()()()()d 0.d b a I y p x y x q x y x f x x x ααηηα'=⎡⎤⎣⎦'+=-+-=⎰因为()x η是任意函数，所以必有()()()()()()0p x y x q x y x f x ''-+-≡. (8.3.8)否则，若在[,]a b 上某点0x 处有()00000()()()()()0p x y x q x y x f x ''-+-≠，不妨设()00000()()()()()0p x y x q x y x f x ''-+->，则由函数的连续性知，在包含0x 的某一区间00[,]a b 上有()()()()()()0p x y x q x y x f x ''-+->.作002200000,[,]\[,],()()(),.x a b a b x x a x b a x b η∈⎧⎪=⎨--≤≤⎪⎩ 显然2()C [,]x a b η∈，且()()0a b ηη==，但()()()()()()()d ba p x y x q x y x f x x x η⎡⎤''-+-⎢⎥⎣⎦⎰ ()00()()()()()()d 0b a p x y x q x y x f x x x η⎡⎤''=-+->⎢⎥⎣⎦⎰，这与式(8.3.7)矛盾. 于是式(8.3.8)成立，即变分问题(8.3.4)的解y 满足微分方程(8.3.1), 且(),()y a y b αβ==故它是边值问题(8.3.1), (8.3.2)的解.(必要性) 设()y y x =是边值问题(8.3.1), (8.3.2)的解，证明y 是变分问题(8.3.4)的解；即：y 使泛函()I y 在21C [,]a b 上达到极小.因为()y y x =满足方程(8.3.1)，所以()()()()()()0p x y x q x y x f x ''-+≡.设任意21()C [,]y y x a b =∈，则函数()()()x y x y x η=-满足条件()()0a b ηη==，且2()C [,]x a b η∈. 于是()()[]()222222()()()()()[()()]()[()()]2()[()()]d ()[()]()[()]2()()d 2()()()()()()()()d ()[()]()[()]d baba baaI y I y I y I y p x y x x q x y x x f x y x x x p x y x q x y x f x y x xp x y x x q x y x x f x x x p x x q x x xηηηηηηηηη-=+-''=+++-+'-+-''=+-++⎰⎰⎰()()()22222()()()()()()d ()[()]()[()]d ()[()]()[()]d .bb ba a bap x y x q x y x f x x x p x x q x x xp x x q x x x ηηηηη⎡⎤'''=--+++⎢⎥⎣⎦'=+⎰⎰⎰⎰因为()0,()0p x q x >≥，所以当()0x η≠时，()22()[()]()[()]d 0bap x x q x x x ηη'+>⎰, 即()()0I y I y ->.只有当()0x η≡时，()()0I y I y -=. 这说明y 使泛函()I y 在21C [,]a b 上达到极小. 证毕！定理8.3.2 边值问题(8.3.1), (8.3.2)存在唯一解.证明用反证法. 若12(),()y x y x 都是边值问题(8.3.1), (8.3.2)的解，且不相等，则由定理8.3.1知，它们都使泛函()I y 在21C [,]a b 上达到极小，因而12()()I y I y > 且 21()()I y I y >，矛盾！因此边值问题(8.3.1), (8.3.2)的解是唯一的.由边值问题解的唯一性，不难推出边值问题(8.3.1), (8.3.2)解的存在性(留给读者自行推导).8.3.2 有限元法等价性定理说明，边值问题(8.3.1), (8.3.2)的解可化为变分问题(8.3.4)的求解问题. 有限元法就是求变分问题近似解的一种有效方法. 下面给出其解题过程：第1步对求解区间进行网格剖分01,i n a x x x x b =<<<<<=区间1[,]i i i I x x -=称为单元，长度1(1,2,,)i i i h x x i n -=-=称为步长，1max i i nh h ≤≤=. 若(1,2,,)i h h i n ==，则称上述网格剖分为均匀剖分.给定剖分后，泛函(8.3.3)可以写成()22()()[()]()()2()()d baI y p x y x q x y x f x y x x '=+-⎰()12211()[()]()()2()()d i i nnx i x i i p x y x q x y x f x y x xS -=='=+-∑∑⎰记. (8.3.9)第2步构造试探函数空间。

第5章两点边值问题求解方法

Page 12
5.2 打靶法作业题5：
用牛顿打靶法求解两点边值问题
y
2 x
y
2 x2
y
sin(ln x) x2 ,
1 x2
y(1) 1, y(2) 2
迭代初始条件取 y(1) 0 。
2020/9/1
航空航天中的计算方法
Page 13
5.3 有限差分法
5.3 有限差分法
用差分近似代替微分，将微
以二阶系统为例，边值问题：分方程化为代数方程求解
y(x) f (x, y(x), y(x)), x a,b
y(a) A, y(b) B
有限差分近似
将区间 xa,b 等分为N个子区间
h ba, N
xi a ih,
i 0,1, 2,
,N
将 y(x)在xi处Taylor展开：
y( x)
y( xi )
航空航天中的计算方法
授课教师：陈琪锋中南大学航空航天学院
第二部分边值问题求解方法
第5章两点边值问题求解方法
内容提要
5.1 常微分方程边值问题的概念 5.2 打靶法 5.3 有限差分法 5.4 有限元法
[1] Part 3: Two-Point Boundary Value Problems.
y( xi )
y( xi1 ) h
y( xi )
一阶精度
若取x=xi-1=x-ih：
y( xi1 )
y( xi )
y( xi )h
1 2
h2 y( xi )
1 6
h3
y(
[2] David L. Darmofal, Computational Methods in Aerospace Engineering (Lecture Notes), MIT, 2005. Chap11,12.

§5.7 边值问题的数值解法

© 2009, Henan Polytechnic University §7 边值问题的数值解法
1010
第五章常微分方程数值解法
即为如下的线性方程组：
1 1 2 h2q 1 1 1 2 h 2 q2 y0 y h2 f 1 1 y2 h2 f 2 ... 1 yn1 h2 f n1 1 yn
© 2009, Henan Polytechnic University §7 边值问题的数值解法
i 1, ... , N 1
6 6
第五章常微分方程数值解法
考虑如下的线性方程
y( x ) p( x ) y( x ) q( x ) y( x ) r ( x ), a x b y(a ) , y(b)
先猜测一个初始斜率 y (a) = s，通过解初值问题 y f ( x , y , y) y(a ) a y(b) = (s) y( a ) s 找出s*使得(s*) = ，即把问题转化为求方程 (s) = 0 的根。
© 2009, Henan Polytechnic University §7 边值问题的数值解法
© 2009, Henan Polytechnic University §7 边值问题的数值解法
1414
考虑如下的二阶常微分方程的边值问题
y( x ) q( x ) y( x ) f ( x ), a x b y(a ) , y(b) 其中q(x)(0)，f(x)在[a，b]上连续，，为常数。
对应的差分问题是：

Matlab 求解化工常微分方程和偏微分方程

y
x
7
8
9
10
x
y(1) Columns 1 through 13
计算值
0 0.0500 0.1000 0.1500 0.2000 0.2500 0.3000 0.3500 0.4000 0.4500 0.5000 0.5500 0.6000 1.0000 1.0526 1.1109 1.1757 1.2479 1.3287 1.4195 1.5218 1.6378 1.7698 1.9210 2.0951 2.2970 Columns 14 through 26 0.6500 0.7000 0.7500 0.8000 0.8500 0.9000 0.9500 1.0000 1.0500 1.1000 1.1500 1.2000 1.2500 2.5329 2.8107 3.1411 3.5381 4.0206 4.6152 5.3598 6.3078 7.5435 9.1928 11.4614 14.7283 19.5991
Results by using ode45()： x y(1) 1.0e+003 * Columns 1 through 13 0 0.0100 0.0200 0.0300 0.0400 0.0500 0.0600 0.0700 0.0800 0.0900 0.1000 0.1100 0.1200 2.0000 0.8788 0.6133 0.4880 0.4181 0.3762 0.3498 0.3328 0.3218 0.3145 0.3097 0.3065 0.3043 Columns 14 through 18 0.1300 0.3029 0.1400 0.3019 0.1500 0.3013 0.1600 0.3009 0.1700 0.3006

二阶常微分方程边值问题的数值解法

摘要本文主要研究二阶常微分方程边值问题的数值解法。

对线性边值问题，我们总结了两类常用的数值方法，即打靶法和有限差分方法，对每种方法都列出了详细的计算步骤和Matlab程序代码，通过具体的算例对这两类方法的优缺点进行了细致的比较。

关键字：常微分方程边值问题；打靶法；差分法；ABSTRACTThis article mainly discusses the numerical methods for solving Second-Order boundary value problems for Ordinary Differential Equations. On the one hand, we review two types of commonly used numerical methods for linear boundary value problems, i.e. shooting method and finite difference method. For each method, we give both the exact calculating steps , we compare the advantages and disadvantages in detail of these two methods through a specific numerical example.Key words：Boundary-Value Problems for Ordinary Differential Equations；Shooting Method；Finite Difference Method；目录第一章引言................................................................................................................... - 1 -第二章二阶线性常微分方程.................................................................................. - 2 -2.1试射法（“打靶”法） ............................................................................................ - 3 -2.1.1简单的试射法............................................................................................ - 3 -2.1.2 基于叠加原理的试射法........................................................................... - 4 -2.2 有限差分法......................................................................................................... - 10 -2.2.1 有限差分逼近的相关概念...................................................................... - 11 -2.2.2 有限差分方程的建立............................................................................. - 13 -2.2.3 其他边值条件的有限差分方程............................................................. - 14 -2.2.4 有限差分方程的解法............................................................................. - 16 -第三章二阶非线性微分方程........................................................ 错误!未定义书签。

常微分方程边值问题有限元

常微分方程边值问题有限元有限元法可归结为如下几个步骤：1.转换成变分问题（应该会用到边界条件）2）对解域进行刨分（可以是不均匀）3）构造基函数（本篇采用基于线性插值的基函数）4）推导出有限元方程5）求解有限元方程6）收敛性和误差估计# 奇怪：这个库必须放在最前面才能一次加载成功usingPlotsgr()# 解决Colab不显示输出数学公式的问题using Markdown: MD, LaTeXfunction latex(expr)expr |> tex |> LaTeXend;考虑两点边值问题：\[ \begin{aligned} &Lu=-\frac{d}{dx}\left(p\frac{du}{dx}\right)+qu=f ,\qquad a < x < b \\ & u(a)=0,\frac{d}{dx}u(b)=0 \end{aligned} \]对应的变分方程(用到两个边界条件)：\[ \begin{aligned}&a(u,v)=(f,v),\qquad \forall v \\ & a(u,v)=\int^b_a{\left(p\frac{du}{dx}\frac{dv}{dx}+quv\ right)dx}\end{aligned} \]区间剖分（可以不均匀）：\[ a = x_0 < x_1 < \dots < x_n=b \\ h_i=x_i-x_{i-1}, \qquad h = \mathrm{max} \{h_i\} \]与此对应的近似解序列 $\{u_i\}$ （待求）为：\[ u_0=0, u_1,u_2,\dots,u_n \]通过这 $n+1$ 个点的值可进行线性插值得到近似解$u_h(x)$ ：\[ \begin{aligned}u_h(x)=\frac{x_i-x}{h_i}u_{i-1}+\frac{x-x_{i-1}}{h_i}u_i,\qquad x_{i-1} < x < x_i \end{aligned} \]可引入线性无关的“山形函数”序列 $\varphi_i(x)$ 构成基底：\[ \boxed{\xi = \frac{x-x_{i-1}}{h_i}\qquad x_{i-1} < x < x_i\qquad 0<\xi<1} \]\[ \varphi_0(x)=\left\{\begin{aligned}&1-\xi, &x_0 < x < x_1 \\ &0, & 其它\end{aligned}\right . \]\[ \varphi_i(x)=\left\{\begin{aligned}&\xi, &x_{i-1} < x < x_i \\ &1-\xi,& x_i < x < x_{i+1} \\ &0, & 其它\end{aligned}\right . \qquad i=1,2,\dots,n-1 \]\[ \varphi_n(x)=\left\{\begin{aligned}&\xi, &x_{n-1} < x < x_n \\ &0, & 其它\end{aligned}\right . \]不难验证，前面的线性插值函数 $u_h(x)$ ，可用这个序列$\varphi_i$ 为基底展开(下面考虑了 $u_0=0$ )：\[ u_h(x)=\sum^n_{i=1}{u_i \varphi_i(x)} \]由图可知（亦可简单验算）：\[ \varphi_i(x)\varphi_j(x)=0, \qquad |i-j|\ge 2 \\\frac{d\varphi_i}{dx}\frac{d\varphi_j}{dx}=0, \qquad|i-j|\ge 2 \]所以 $a(\varphi_i,\varphi_j)$ 的非零值只可能有：\[ \begin{aligned}&a(\varphi_{1},\varphi_1),a(\varphi_ {2},\varphi_1) \\ &a(\varphi_{j-1},\varphi_j),a(\varphi_{j},\varphi_j),a(\varphi_{j+1} ,\varphi_j),\qquad j=2,\dots,n-1\\ & a(\varphi_{n-1},\varphi_n),a(\varphi_{n},\varphi_{n}) \end{aligned} \]可分别算出这三类非零值：\[ a(\varphi_{j-1},\varphi_j)=\int^1_0{\left[-h_j^{-1}p(x_{j-1}+h_j\xi)+h_j q(x_{j-1}+h_j\xi)(1-\xi)\xi\right]d\xi} \]\[ \begin{aligned}a(\varphi_j,\varphi_j)&=\int^1_0{\le ft[h_j^{-1}p(x_{j-1}+h_j\xi)+h_j q(x_{j-1}+h_j\xi)\xi^2\right]d\xi} \\ &+\int^1_0{\left[h_{j+1}^{-1}p(x_{j}+h_{j+1}\xi)+h_{j+1} q(x_{j} +h_{j+1}\xi)(1-\xi)^2\right]d\xi}\end{aligned} \] \[ a(\varphi_{j+1},\varphi_j)=\int^1_0{\left[-h_{j+1}^{-1}p(x_{j}+h_{j+1}\xi)+h_{j+1}q(x_{j}+h_{j+1}\xi)(1-\xi)\xi\right]d\xi} \]\[ a(\varphi_n,\varphi_n)=\int^1_0{\left[h_n^{-1}p(x_{n-1}+h_n\xi)+h_n q(x_{n-1}+h_n\xi)\xi^2\right]d\xi} \]using SciPy# # “山形”基函数#functio n φ(t,i,x,h,n)if i >1&& t>=x[i-1] && t<=x[i]return (t-x[i-1])/h[i-1]elseif i < n+2&& t>x[i] && t<=x[i+1]return1-(t-x[i])/h[i]elsereturn0endend;# # a(φ_i,φ_j)## 假设已有p(t),q(t)#functionaa(i,j,x,h,n)if abs(i-j) >=2return0;end;if i==j-1α(ξ)=-p(x[j-1]+h[j-1]*ξ)/h[j-1]+h[j-1]*q(x[j-1]+h[j-1]*ξ)*(1-ξ)*ξ;return SciPy.integrate.quad(α, 0, 1)[1];end;if i==j && j<n+1β(ξ)=p(x[j-1]+h[j-1]*ξ)/h[j-1]+h[j-1]*q(x[j-1]+h[j-1]*ξ)*ξ^2+p(x[j]+h[j]*ξ)/h[j]+h[j]*q(x[j]+h[j]*ξ)* (1-ξ)^2;return SciPy.integrate.quad(β, 0, 1)[1];end;if i==j && j==n+1βn(ξ)=p(x[j-1]+h[j-1]*ξ)/h[j-1]+h[j-1]*q(x[j-1]+h[j-1]*ξ)*ξ^2;return SciPy.integrate.quad(βn, 0, 1)[1];end;if i==j+1γ(ξ)=-p(x[j]+h[j]*ξ)/h[j]+h[j]*q(x[j]+h[j]*ξ)*(1-ξ)*ξ;return SciPy.integrate.quad(γ, 0, 1)[1];end;end;正式求解的julia代码：using Random# 参数选定a=0; b=1f(t)=(pi^2/2)sin(pi*t/2);p(t) =1;q(t) = pi^2/4;function NSolve(n)# 随机生成的不均匀的区间刨分h = [0.5*rand()+0.5for i in1:n];h = (b-a)*h/sum(h);x = zeros(n+1);x[1]=a; x[n+1]=b;x[2:n]=[sum(h[1:i]) for i in1:(n-1)];h[n]=x[n+1]-x[n];# 待求近似解u = zeros(n+1);# 计算a(φ_i,φ_j)矩阵A = [aa(i,j,x,h,n) for j in2:n+1, i in2:n+1];# 计算(f,φ_j)列向量fφ = [SciPy.integrate.quad(t->f(t)*φ(t,j,x,h,n), x[j-1],(j==n+1? b : x[j+1]))[1] for j in2:n+1];# 解出插值端点值u[2:n+1] = in v(A)*fφ;# 近似解序列(x,u)end;做为对比，也可解出精确解：using SymPy@vars ty = SymFunction("y")diffeq =-diff(y(t),t,2)+(PI^2/4)y(t)⩵(PI^2/2)*sin(t*PI/2)# 通解ex = dsolve(diffeq, y(t))# 根据边界条件确定积分常数ex1 = rhs(ex)ex2 = diff(ex1,t)eqs = [ex1(t=>0)⩵0,ex2(t=>1)⩵0]C = solve(eqs,[Sym("C1"),Sym("C2")])# 解析解ex = ex(C)ex |> latex |> MD\[ y{\left(t \right)} = \sin{\left(\frac{\pi t}{2} \right)} \]绘制近似解和精确解# 精确解plot(rhs(ex),color="black",xlims=(-0.1,1.1),label="U",xlabel="x",ylabel="u(x)")# 近似解 n=4,8plot!(NSolve(4),color="green",label="u4")plot!(NSolve(8),color="red",label="u8")《微分方程数值解法（第4版）-李荣华＆刘播-高等教育出版社-2009》。

常系数微分方程组的解法

幂级数法
将高阶线性微分方程转化为幂级数形式，然后通过幂级数的性质求解方程。
高阶非线性微分方程的解法
分离变量法
将非线性微分方程转化为多个一阶微分方程，然后分别求解。
迭代法
通过迭代公式逐步逼近非线性微分方程的解。
数值解法
利用数值计算方法求解非线性微分方程的近似解，如欧拉法、龙格-库塔法等。
05
解决微分方程组对于理解复杂系统的行为和预测未来发展趋势具有重要意义。
常系数微分方程组的定义
常系数微分方程组是指方程中的系数为常数的一类微分方程组。
常系数微分方程组的一般形式为 dy/dx = f(x, y)，其中 f(x, y) 是已知的函数。
02
线性常系数微分方程组的解法
特征根法
总结词
神经传导
在神经传导过程中，微分方程组可以用来描述神经信号的传递速度和传导通路的建立。
生态系统的稳定性
微分方程组可以用来分析生态系统的稳定性，如物种之间的相互作用和生态平衡的维持。
THANKS
感谢观看
特征根法是一种通过解方程的特征方程来求解线性常系数微分方程组的方法。
详细描述
特征根法的基本思想是，对于形如$y'' + py' + qy = 0$的一阶线性常系数微分方程，通过求解其特征方程$lambda^2 + plambda + q = 0$，得到其特征根$lambda_1$和 $lambda_2$，然后利用这些特征根来求解原微分方程。
线性微分方程的方法。
02
通过将多个变量分离，可以将一个复杂的微分方程组
分解为多个简单的微分方程，从而简化求解过程。
03

关于常微分方程组边值问题的研究的开题报告

关于常微分方程组边值问题的研究的开题报告尊敬的导师：您好！感谢您在百忙之中抽出时间阅读我的开题报告，本文将简述我对常微分方程组边值问题的研究的初步规划。

1.研究的背景和意义常微分方程组是数学的重要分支之一，边值问题是研究常微分方程组的重要问题之一。

常微分方程组边值问题广泛应用于物理、化学、力学等各个领域，在实际应用中有非常重要的作用。

因此，对于常微分方程组边值问题的深入研究对于学术界和工业界都具有重要的意义。

2.研究的目的和内容本次研究的目的在于进一步深入研究常微分方程组边值问题，探究其基本理论。

我们将从以下方面进行研究：（1）研究边值问题的定义和性质。

（2）研究常微分方程组的解的存在性和唯一性。

（3）研究解的连续依赖性和解的依赖于参数的连续性。

（4）研究常微分方程组边值问题的边值条件。

3.研究方法和技术路线我们将采用数学分析的方法，结合常微分方程理论和边值问题的基本概念进行研究。

具体的技术路线如下：（1）建立常微分方程组的数学模型。

（2）引入边值问题的相关概念，对问题进行分析。

（3）探究解的存在性和唯一性。

（4）研究解的连续依赖性和解的依赖于参数的连续性。

（5）研究边界条件。

（6）进行数学分析并给出定理和解的表达式。

4.预计成果我们将通过以上的研究方法和技术路线，得出如下成果：（1）较为系统和深入地掌握常微分方程组边值问题的基本理论和方法。

（2）得出常微分方程组边值问题的一些重要性质和解的表达式。

（3）通过分析实际问题，得出一些有用的结论和应用方法。

5.研究计划本次研究预计耗时一年，具体的研究计划如下：第一季度：研究边值问题的基本概念和定义；第二季度：研究常微分方程组的解的存在性和唯一性；第三季度：研究解的连续依赖性和解的依赖于参数的连续性；第四季度：研究常微分方程组边值问题的边值条件。

6.结语本次研究将会进一步深入研究常微分方程组边值问题的基本理论和方法，得出一些有用的结论和应用方法。

希望在导师的指导下能够顺利完成研究任务，有益于学术界和工业界。