常微分方程组(边值)
重要：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。
常微分方程的边值问题
常微分方程的边值问题一、引言在数学中，微分方程是研究自然界中变化和发展的重要工具。
它描述了物体在不同变化条件下的行为规律，并被广泛应用于物理、工程、经济等领域。
边值问题是微分方程中的一个重要分支，它关注的是在一定边界条件下的解。
二、常微分方程常微分方程是指只含有关于一个自变量的一阶或高阶导数的方程。
一般形式为：[F(x, y, y’, y’’, , y^{(n)}) = 0]其中，x是自变量，y是未知函数。
常微分方程的求解可以分为两种类型：初值问题和边值问题。
三、边值问题的定义边值问题是指在一定边界条件下，求解微分方程的解。
对于二阶常微分方程，边值问题的一般形式为：[y’‘(x) = f(x, y, y’)， a < x < b， y(a) = , y(b) = ]其中，a和b是给定的边界点，()和()是给定的边界值。
四、边值问题的求解方法边值问题的求解可以分为两种方法：迭代方法和直接方法。
4.1 迭代方法迭代方法是通过不断迭代逼近的方式求解边值问题。
常用的迭代方法有有限差分法和有限元法。
4.1.1 有限差分法有限差分法是一种将微分方程转化为差分方程进行求解的方法。
它将求解域离散化，并通过差分近似来近似微分项，最终通过迭代逼近求得边界值。
有限差分法的基本思想是将求解域划分为若干个离散的网格点，然后使用近似公式将微分项替换为差分项，从而得到差分方程。
通过迭代求解差分方程，最终得到边界条件下的解。
4.1.2 有限元法有限元法是一种将微分方程转化为代数方程组进行求解的方法。
它通过将求解域划分为有限个小区域，然后在每个小区域上选择一个试验函数来代表解，在满足边界条件的情况下，通过最小化误差的方法得到近似解。
有限元法的基本思想是将求解域划分为若干个小单元，然后在每个小单元上选择一个适当的试验函数，通过建立弱形式和加权残差方法得到代数方程组，最终通过迭代求解代数方程组得到边界条件下的解。
4.2 直接方法直接方法是通过对微分方程进行直接求解的方法，其中最常用的方法是变分法。
常微分方程的边值问题
常微分方程的边值问题常微分方程是数学中一个重要的分支，研究的是函数的导数与自变量之间的关系。
在实际问题中，常微分方程的解可以描述物理、工程、经济等领域的变化规律。
而边值问题是常微分方程中的一类特殊问题，它要求在给定的边界条件下求解方程的解。
一、边值问题的定义与分类边值问题是指在一定边界条件下求解常微分方程的解。
边界条件是一组给定的条件，它们通常是关于未知函数及其导数在一些特定点上的值或关系。
边值问题可分为以下两类：1. Dirichlet 边值问题：给定函数在边界上的值。
假设我们要求解的常微分方程为 y''(x) + p(x)y'(x) + q(x)y(x) = r(x)，边值问题可以表示为：y(a) = A，y(b) = B其中，a, b 是给定的自变量取值，A, B 是给定的常数。
2. Neumann 边值问题：给定函数在边界上的导数值。
假设我们要求解的常微分方程还是 y''(x) + p(x)y'(x) + q(x)y(x) = r(x)，边值问题可以表示为：y'(a) = A，y'(b) = B二、求解边值问题的方法求解边值问题有多种方法，其中比较常用的包括：1. 分离变量法这是一种基本的求解边值问题的方法。
通过将方程中的未知函数分离变量，得到一个关于自变量的方程和一个关于未知函数的方程，再分别求解这两个方程。
2. 特征值法对于某些特殊的边值问题，可以使用特征值法进行求解。
特征值法的关键在于将边值问题转化为一个特征值问题，通过求解特征值和特征函数来得到方程的解。
3. 迭代法对于某些复杂的边值问题，可以使用迭代法逐步逼近方程的解。
迭代法是通过不断逼近函数解来改善近似解的精度，从而得到较为准确的解。
三、常见的边值问题应用常微分方程的边值问题在实际应用中具有广泛的应用，下面列举几个常见的例子：1. 自由振动问题自由振动是常微分方程的一个典型应用，比如弹簧振子的运动可以用一阶线性常微分方程来描述。
常微分方程边值问题的数值解法
……
……
……
3.6896236 3.6865656 -3.058*10-3
4.5635316 4.5612288 -2.303*10-3
5.5854269 5.5841425 -1.284*10-3
6.7725887 6.7725887 0
数值计算方法
对区间[a,b]作等距分划: x j a jh( j 0,1,2,...n)
h b a。由数值微分公式 n
y(x j )
y(x j1) y(x j1) 2h
h2 6
y(3) ( j )
y(x j )
y(x j1) 2 y(x j ) h2
y(x j1)
h2 12
y(4) ( j )
其中, j , j (x j1, x j1)。
差分方程的建立
代入y p(x) y q(x) y r(x), x [a,b]得差分方程：
y j1 2 y j h2
y j1
pj
y j1 y j1 2h
qj yj
rj ( j 1,2,...n 1)
这是求y j ( j 0,1,2,...,n)的n 1个方程,还缺的两个方
程由边界条件补出。
差分方程的建立
对于第一类边界条件：y0 , yn ,即已给出
两个未知量的解, 这时整理后有
b1 c1 a2 b2 c2 y1 d1 a1 源自y2d2an2
bn2
cn2
yn2
dn2
an1 bn1 yn1 dn1 cn1
其中a
解 : 此方程的解析解为y x2 x3 1 x4 x2 ln x。 2
例题
现取步长h 0.1，此时p(x) 2 ， x
课件：级第四章 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的初值问题为-版凤防行：小人学数：学烧
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
常微分方程组边值问题解法打靶法Shooti ng Method (shoot in g.m )% 丁靶法求常微分方程的边值问题function [x,a,b ,n]=shooti ng(fu n, xO,x n, eps) if nargin<3eps=1e-3;endx1=x0+ra nd;[a,b]=ode45(fu n, [0,10],[0,x0]');c0=b(le ngth(b),1);[a,b]=ode45(fu n, [0,10],[0,x1]');c1=b(le ngth(b),1);x2=x1-(c1-x n)*(x1-x0)/(c1-c0);n=1;while (no rm(c1-x n)>=eps & no rm(x2-x1)>=eps) x0=x1;x 仁x2;[a,b]=ode45(fu n,[ 0,10],[0,x0]');cO=b(le ngth(b),1);[a,b]=ode45(fu n,[ 0,10],[0,x1]');c1= b(le ngth(b),1)x2=x1-(c1-x n)*(x1-x0)/(c1-c0);n=n+1;endx=x2;应用打靶法求解下列边值问题:y 10 0 解：将其转化为常微分方程组的初值问题命令:xO=[O:O.1:1O];y0=32*((cos(5)-1)/si n( 5)*si n(x0/2)-cos(x0/2)+1); plot(xO,yO,'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))dy i dx y 2 dy 2 dx y i 0y 4 y odydx X0真实解30'12^4567^9 10函数：(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]=shoot in g('odebvp',10,0,1e-3)计算结果：(eps=0.001 )t=11.9524plot(x,y(:,1))x0=[0:1:10];y0=32*((cos(5)-1)/si n( 5)*si n(x0/2)-cos(x0/2)+1); hold onplot(xO,yO, ' o')有限差分法 Finite Differenee Methods FDM同上例:Y i i2Y i y i ih 2 2 ——14函数：(differe nce.m )%有限差分法求常微分方程的边值问题 function [x,y]=differe nce(xO,x n,yO,yn,n) h=(x n-xO)/n;a=eye( n-1)*(-(2-h A 2/4)); for i=1: n-2 a(i,i+1)=1; a(i+1,i)=1; endb=o nes( n-1,1)*8*hA2; b(1)=b(1)-0; b(n-1)=b( n-1)-0; yy=a\b;x(1)=x0;y(1)=y0; for i=2: n x(i)=x0+(i-1)*h;Y i i 2若划分为10个区间，则：h- Y i Y i 1 8h 24(differe nce.m )1h 2 2 411Y 1 Y 28h 28h 20 h 2Y n 2 8h 224 1Y n 18h 2 0.2h 124y(i)=yy(i-i);endx(n)=xn;y(n)=yn;命令：[x,y]=differe nce(0,10,0,0,100);计算结果：xO=[O:O.1:1O];y0=32*((cos(5)-1)/si n( 5)*si n(x0/2)-cos(x0/2)+1);plot(xO,yO,'r')hold on[x,y]=differe nce(0,10,0,0,5);plot(x,y,'.')[x,y]=differe nce(0,10,0,0,10);plot(x,y,'--')[x,y]=differe nce(0,10,0,0,50);plot(x,y,'-.')正交配置法Orthogonal Collocatioin Methods CM构造正交矩阵函数(collmatrix.m )%E交配置矩阵(均用矩阵法求对称性与非对称性正交配置矩阵) function [am,bm,wm,an,bn,wn]=collmatrix(a,m,fm,n,fn) 真实解的3 4 6 7 9 9 1DXxO=symm(a,m,fm); %a为形状因子;m为零点数;fm为对称的权函数(0为权函数1,非0为权函数1-x A2)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)A(2*i-2);cm(j,i)=(2*i-2)*xm(j)A(2*i-3);dm(j,i)=(2*i-2)*(2*i-3+(a-1))*xm(j)A(2*i-3+(a-1)-1-(a-1));endfmm(j)=1/(2*j-2+a);endam=cm*i nv(qm);bm=dm*i nv(qm);wm=fmm*i nv(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)=x n(j)A(i-1);if j==0 | i==1cn (j,i)=0;elsecn (j,i)=(i-1)*x n(j)A(i-2);endif j==0 | i==1 | i==2dn (j,i)=O;elsedn (j,i)=(i-2)*(i-1)*x n(jF(i-3);endendfnn (j)=1/j;endan=cn *i nv(qn);bn=dn *i nv(qn); wn=fnn*inv(qn); %E交多项式求根(适用于对称问题)fun ctio n p=symm(a,m,fm) %a 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);% elsec2=c2*(m+a/2+j+1);% endc3=c3*(a/2+j);c4=c4*(1+j);endp(m+1-i)=c1*c2/c4/c3;endp(m+1)=1;%为多项式系数向量p=sqrt(roots(p)); 为形状因子,m为配置点数,fm为权函数权函数为1权函数为1-x A2,求出根后对对称问题还应开方才是零点%E交多项式求根（适用于非对称性问题）解：（1）标准化令x r/R ，y C/C S 代入微分方程及边界条件得：丄2 x 2dy36yx 2dx dxx 0単0 dx x 1, y 1fun ctio n p=un sym m(n,fn) if fn==0r(1)=(-1)A n*n *( n+1);% 权函数为 1elser(1)=(-1)A n*n*(n+2);% 权函数为 1-x end for i=1: n-1 if fn==0r(i+1)=( n-i)*(i+n+1)*r(i)/(i+1)/(i+1);% elser(i+1)=( n-i)*(i+n+2)*r(i)/(i+1)/(i+1);% end end for j=1: np( n+1-j)=(-1)A(j+1)*r(j); endp(n +1)=(-1)A( n+1); p=roots(p);权函数为1权函数为1-x 应用正交配置法求解以下等温球形催化剂颗粒内反应物浓度分布，其浓度分布的数学 模型为：丄r 2 drr r2dC r dr36C R 2dC 0，丁 dr(2)离散化N 1B ji y i 36y j 0i 1j 1, 2,N 1(3)转化为代数方程组(以 N 3为例)因为y N i 目41，所以整理上式得:B 11 36 B 12 B 13B 21 B 22 36 B 23 B 31 B 32 B 33 36B 41 B 42 B 43y iy 2 y 3 B 14 B 24B 34 B 44 36 本例中的代数方程组为线性方程组，可采用线性方程组的求解方法；若为非线性方程组则 采用相应的方法求解。
命令：N=3,权函数为1-x [am,bm,wm,a n,bn,wn ]=collmatrix(3,3,1,3,1);b1=bm;for i=1:4 b1(i,i)=bm(i,i)-36; end a0=b1(1:4,1:3); b0=-b1(1:4,4); y=aO\bO; y(4)=1;p=exam31(3,3);(注意要对文件修改权函数为 x=[0.3631,0.6772,0.8998,1]; %零点 plot(x,y,'o')hold onx0=0:0.1:1; % 真实解 y0=si nh(6*x0)./x0/si nh(6); plot(xO,yO,'r') (只用对称性配置矩阵) 1-x 2)若权函数改为1,则以下语句修改，其他不变 [am,bm,wm,a n,b n,w n]=collmatrix(3,3,0,3,1);(只用对称性配置矩阵)r1 2 3 4 yyyy3^142434BBB33433421112 3 4BBBp=exam31(3,3);(注意要对文件修改权函数为1)x=[0.4058,0.7415,0.9491,1]; % 零点计算结果:权函数为1- x权函数为1%11正交配置法0.9真实解0.80.7 ・0.6 -y 0.5 -0.4 - -0.3 -X：0.2 --0.1 --1ii 1 ■0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1x边值问题的MatLab解法精确解:2x函数：(collfuni.m ) fun ctio n f=collfu ni(x,y) f(i)=y(2);f(2)=4*y(1);f=[f(1);f(2)];(collbci.m ) function f=collbc1(a,b) f=[a(1)-1;b(1)-exp(2)];命令：soli ni t=bvpi ni t([0:0.1:1],[1,1]) sol=bvp4c(@collfu n1,@collbc1,soli nit) plot(sol.x,sol.y)hold onplot(sol.x,exp(2*sol.x),'*')yi y2 真实解y 4y 0 x 1y 0 1, y 1 e2y i yy2 4y iy i 0 i, y i i e2真实解精确解:函数：(collfun2.m )function f=collfun2(x,y) f(1)=y(2)；f(2)=(1-x.A 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)];命令：soli nit=bvpi ni t([0:0.1:1],[1,1]);sol=bvp4c(@collfu n2,@collbc2,soli nit); plot(sol.x,sol.y) hold onplot(sol.x,(sol.x-1).*exp(-sol.x),'*')y x 1 y 2y 1 x 2 e x0 x 1y 0 2, y 1 1/ey i y 2y 2 1 x 2 e x 2y 1 % 0 2, y 1 1x 1 y 21/e 真实解函数：(collfu n3.m ) function f=collfun3(x,y) f(i)=y(2);f(2)=(2-log(x))./x+y(1)./x-y(2).A2; f=[f(i)；f(2)]；(collbc3.m )function f=collbc3(a,b) f=[2*a(i)-a (2);b(2)-1.5];命令：soli nit=bvpi ni t([1:0.1:2],[1,1]);sol=bvp4c(@collfu n3,@collbc3,soli nit); plot(sol.x,sol.y) hold onplot(sol.x,sol.x+log(sol.x),'*')真实解2y yy12ln x 1 x 2x x2y 1y 1 0,y 23/2精确解：y x In xy 2 2 In xx2y i 1 y 2 i 0,y i y 2xy 2 23/2在260 C 的基础面上，为促进传热在此表面上 增加纯铝的圆柱形肋片，其直径为25mm 高为150mm 该柱表面受到16 C 气流的冷却，气流与肋 片表面的对流传热系数为15W/m 2 K ，肋端绝热;肋片的导热系数为 236W/m K ，假设肋片的导热 热阻与肋片表面的对流传热热阻相比可以忽略；试 求肋片中的温度分布，及单个肋片的散热量为多 少？解：根据以上条件可知：导热热阻与对流热阻相比可以忽略，则在肋片径向上没有温 度分布，在轴向上存在温度分布，外界气流与肋片的对流传热则可转化为内热源；故该问 题为导热系数为常数的一维稳定热传导，其导热微分方程为：这是两点边值的常微分方程求解问题，故可转化为如下形式:y iy 2hp y yA Cy 2 H 0, 0边界条件为：x 0时，t 0 分析解：t t t 0 td 2t dx 2hp t t A C260 C （肋根）； x H 时，dtdx0 （肋端绝热）。