常微分方程组(边值)
重要：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转向③继续计算直到满意为止。
特别地，若微分方程是线性的，则打靶法变成线性组合法，即根据常微分方程理论适当选取初值可得到一组线性独立解，利用它们的线性组合导出边值问题的解。
例如线性方程边值问题的数值解可通过两个初值问题数值解来实现。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
常微分方程组边值问题解法打靶法Shooting Method （shooting.m ）%打靶法求常微分方程的边值问题function [x,a,b,n]=shooting(fun,x0,xn,eps)if nargin<3eps=1e-3;endx1=x0+rand;[a,b]=ode45(fun,[0,10],[0,x0]');c0=b(length(b),1);[a,b]=ode45(fun,[0,10],[0,x1]');c1=b(length(b),1);x2=x1-(c1-xn)*(x1-x0)/(c1-c0);n=1;while (norm(c1-xn)>=eps & norm(x2-x1)>=eps)x0=x1;x1=x2;[a,b]=ode45(fun,[0,10],[0,x0]');c0=b(length(b),1);[a,b]=ode45(fun,[0,10],[0,x1]');c1=b(length(b),1)x2=x1-(c1-xn)*(x1-x0)/(c1-c0);n=n+1;endx=x2;应用打靶法求解下列边值问题：()()⎪⎪⎩⎪⎪⎨⎧==-=010004822y y y dxy d解：将其转化为常微分方程组的初值问题()⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧==-==t dx dy y y y dx dy y dx dy x 0011221048命令：x0=[0:0.1:10];y0=32*((cos(5)-1)/sin(5)*sin(x0/2)-cos(x0/2)+1); 真实解plot(x0,y0,'r')hold on[x,y]=ode45('odebvp',[0,10],[0,2]');plot(x,y(:,1))[x,y]=ode45('odebvp',[0,10],[0,5]');plot(x,y(:,1))[x,y]=ode45('odebvp',[0,10],[0,8]');plot(x,y(:,1))[x,y]=ode45('odebvp',[0,10],[0,10]');plot(x,y(:,1))函数：（odebvp.m）%边值常微分方程(组)函数function f=odebvp(x,y)f(1)=y(2);f(2)=8-y(1)/4;f=[f(1);f(2)];命令：[t,x,y,n]=shooting('odebvp',10,0,1e-3)计算结果：（eps=0.001）t=11.9524plot(x,y(:,1))x0=[0:1:10];y0=32*((cos(5)-1)/sin(5)*sin(x0/2)-cos(x0/2)+1); hold onplot(x0,y0,’o’)有限差分法Finite Difference Methods FDM （difference.m ）同上例：4822y dx y d -=⇒482211i i i i y h y y y -=+--+ 2121842h y y h y i i i =+⎪⎪⎭⎫ ⎝⎛--+- 若划分为10个区间，则： ⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡--=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⎪⎪⎭⎫ ⎝⎛--⎪⎪⎭⎫ ⎝⎛--⎪⎪⎭⎫ ⎝⎛--⎪⎪⎭⎫ ⎝⎛----08880842114211421142222212212222h h h h y y y y h h h h n n M M O O O函数：（difference.m ）%有限差分法求常微分方程的边值问题function [x,y]=difference(x0,xn,y0,yn,n)h=(xn-x0)/n;a=eye(n-1)*(-(2-h^2/4));for i=1:n-2a(i,i+1)=1;a(i+1,i)=1;endb=ones(n-1,1)*8*h^2;b(1)=b(1)-0;b(n-1)=b(n-1)-0;yy=a\b;x(1)=x0;y(1)=y0;for i=2:nx(i)=x0+(i-1)*h;y(i)=yy(i-1);endx(n)=xn;y(n)=yn;命令：[x,y]=difference(0,10,0,0,100);计算结果：x0=[0:0.1:10];y0=32*((cos(5)-1)/sin(5)*sin(x0/2)-cos(x0/2)+1); 真实解plot(x0,y0,'r')hold on[x,y]=difference(0,10,0,0,5);plot(x,y,’.’)[x,y]=difference(0,10,0,0,10);plot(x,y,’--’)[x,y]=difference(0,10,0,0,50);plot(x,y,’-.’)正交配置法Orthogonal Collocatioin Methods CM构造正交矩阵函数（collmatrix.m）%正交配置矩阵(均用矩阵法求对称性与非对称性正交配置矩阵)function [am,bm,wm,an,bn,wn]=collmatrix(a,m,fm,n,fn)x0=symm(a,m,fm); %a为形状因子;m为零点数;fm为对称的权函数(0为权函数1,非0为权函数1-x^2)for i=1:mxm(i)=x0(m+1-i);endxm(m+1)=1;for j=1:m+1for i=1:m+1qm(j,i)=xm(j)^(2*i-2);cm(j,i)=(2*i-2)*xm(j)^(2*i-3);dm(j,i)=(2*i-2)*(2*i-3+(a-1))*xm(j)^(2*i-3+(a-1)-1-(a-1));endfmm(j)=1/(2*j-2+a);endam=cm*inv(qm);bm=dm*inv(qm);wm=fmm*inv(qm);x1=unsymm(n,fn); %n为零点数;fn为非对称的权函数(0为权函数1,非0为权函数1-x) xn(1)=0;for i=2:n+1xn(i)=x1(n+2-i);endxn(n+2)=1;for j=1:n+2for i=1:n+2qn(j,i)=xn(j)^(i-1);if j==0 | i==1cn(j,i)=0;elsecn(j,i)=(i-1)*xn(j)^(i-2);endif j==0 | i==1 | i==2dn(j,i)=0;elsedn(j,i)=(i-2)*(i-1)*xn(j)^(i-3);endendfnn(j)=1/j;endan=cn*inv(qn);bn=dn*inv(qn);wn=fnn*inv(qn);%正交多项式求根(适用于对称问题)function p=symm(a,m,fm) %a为形状因子,m为配置点数,fm为权函数for i=1:mc1=1;c2=1;c3=1;c4=1;for j=0:i-1c1=c1*(-m+j);if fm==0c2=c2*(m+a/2+j);%权函数为1elsec2=c2*(m+a/2+j+1);%权函数为1-x^2endc3=c3*(a/2+j);c4=c4*(1+j);endp(m+1-i)=c1*c2/c4/c3;endp(m+1)=1;%为多项式系数向量,求出根后对对称问题还应开方才是零点p=sqrt(roots(p));%正交多项式求根(适用于非对称性问题)function p=unsymm(n,fn)if fn==0r(1)=(-1)^n*n*(n+1);%权函数为1elser(1)=(-1)^n*n*(n+2);%权函数为1-xendfor i=1:n-1if fn==0r(i+1)=(n-i)*(i+n+1)*r(i)/(i+1)/(i+1);%权函数为1elser(i+1)=(n-i)*(i+n+2)*r(i)/(i+1)/(i+1);%权函数为1-xendendfor j=1:np(n+1-j)=(-1)^(j+1)*r(j);endp(n+1)=(-1)^(n+1);p=roots(p);应用正交配置法求解以下等温球形催化剂颗粒内反应物浓度分布，其浓度分布的数学模型为：⎪⎪⎪⎩⎪⎪⎪⎨⎧=====⎪⎭⎫ ⎝⎛S C C r dr dC r R C dr dC r dr d r ,10,0361222 解：（1）标准化令R r x /=，S C C y /=代入微分方程及边界条件得：⎪⎪⎪⎩⎪⎪⎪⎨⎧=====⎪⎭⎫ ⎝⎛1,10,036122y x dx dy x y dx dy x dx d x（2）离散化03611=-∑+=j N i i ji y y B1,2,1+=N j Λ（3）转化为代数方程组（以3=N 为例）⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡----000036363636432144434241343332312423222114131211y y y y B B B B B B B B B B B B B B B B 因为141==+y y N ，所以整理上式得：⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡+----=⎥⎥⎦⎤⎢⎢⎣⎡⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡---3636363644342414321434241333231232221131211B B B B y y y B B B B B B B B B B B B 本例中的代数方程组为线性方程组，可采用线性方程组的求解方法；若为非线性方程组则采用相应的方法求解。
命令：N=3,权函数为1-x2 [am,bm,wm,an,bn,wn]=collmatrix(3,3,1,3,1);（只用对称性配置矩阵）b1=bm;for i=1:4b1(i,i)=bm(i,i)-36;enda0=b1(1:4,1:3);b0=-b1(1:4,4);y=a0\b0;y(4)=1;p=exam31(3,3);（注意要对文件修改权函数为1-x 2）x=[0.3631,0.6772,0.8998,1]; %零点plot(x,y,'o')hold onx0=0:0.1:1; %真实解y0=sinh(6*x0)./x0/sinh(6);plot(x0,y0,'r')若权函数改为1，则以下语句修改，其他不变[am,bm,wm,an,bn,wn]=collmatrix(3,3,0,3,1);（只用对称性配置矩阵）p=exam31(3,3);（注意要对文件修改权函数为1）x=[0.4058,0.7415,0.9491,1]; %零点计算结果：权函数为1- x2权函数为1y 正交配置法真实解边值问题的MatLab 解法()()⎩⎨⎧==<<=''21,10104e y y x y y ⇒ ()()⎪⎩⎪⎨⎧==='='21112211,104ey y y y y y 精确解：x e y 2=函数：（collfun1.m ）function f=collfun1(x,y) f(1)=y(2); f(2)=4*y(1); f=[f(1);f(2)]; （collbc1.m ）function f=collbc1(a,b) f=[a(1)-1;b(1)-exp(2)];命令：solinit=bvpinit([0:0.1:1],[1,1])sol=bvp4c(@collfun1,@collbc1,solinit) plot(sol.x,sol.y) hold onplot(sol.x,exp(2*sol.x),'*') 真实解()()()()⎩⎨⎧='='<<-=-'++''-e y y x e x y y x y x /11,20101212 ⇒ ()()()()⎪⎩⎪⎨⎧='='+-+-='='-e y y y x y e x y y y x /11,2012111212221 精确解：()x e x y --=1函数：（collfun2.m ）function f=collfun2(x,y) f(1)=y(2);f(2)=(1-x.^2).*exp(-x)+2*y(1)-(x+1).*y(2); f=[f(1);f(2)]; （collbc2.m ）function f=collbc2(a,b) f=[a(2)-2;b(2)-exp(-1)];命令：solinit=bvpinit([0:0.1:1],[1,1]);sol=bvp4c(@collfun2,@collbc2,solinit); plot(sol.x,sol.y) hold onplot(sol.x,(sol.x-1).*exp(-sol.x),'*') 真实解()()()()⎪⎩⎪⎨⎧='='-<<-=-'+''2/32,011221ln 212y y y x x xx y y y ⇒ ()()()()⎪⎪⎩⎪⎪⎨⎧==--+-='='2/32,0112ln 21221221221y y y y x y x x y y y 精确解：x x y ln +=函数：（collfun3.m ）function f=collfun3(x,y) f(1)=y(2);f(2)=(2-log(x))./x+y(1)./x-y(2).^2; f=[f(1);f(2)]; （collbc3.m ）function f=collbc3(a,b) f=[2*a(1)-a(2);b(2)-1.5];命令：solinit=bvpinit([1:0.1:2],[1,1]);sol=bvp4c(@collfun3,@collbc3,solinit); plot(sol.x,sol.y) hold onplot(sol.x,sol.x+log(sol.x),'*') 真实解在260C ︒的基础面上，为促进传热在此表面上增加纯铝的圆柱形肋片，其直径为25mm ，高为150mm ；该柱表面受到16C ︒气流的冷却，气流与肋片表面的对流传热系数为15K m W ⋅2/，肋端绝热；肋片的导热系数为236K m W ⋅/，假设肋片的导热热阻与肋片表面的对流传热热阻相比可以忽略；试求肋片中的温度分布，及单个肋片的散热量为多少？解：根据以上条件可知：导热热阻与对流热阻相比可以忽略，则在肋片径向上没有温度分布，在轴向上存在温度分布，外界气流与肋片的对流传热则可转化为内热源；故该问题为导热系数为常数的一维稳定热传导，其导热微分方程为：()C A t t hp dxt d λλ∞-=Φ=&22 边界条件为：0=x 时，C 2600︒=t （肋根）；H x =时，0==Hx dxdt（肋端绝热）。