解各种方程.txt

matlab中解方程的主要命令,和详细的使用方法,希望对朋友门有用!

  > > crit=xd(find(product<0)) % 注意用到find指令 

若上述多项式以中央差分方式计算其函数微分项，就不能以diff计算，而须自行计算如 

下 
： 

  > > num=f(3:length(f))-f(1:length(f)-2); % 注意中央差分是 f(k+1)-f(k-1) 
  > > deno=x(3:length(f))-x(1:length(f)-2); % 注意中央差分是 x(k+1)-x(k-1) 
  > > df_c=num./deno; 
  > > xd=x(2:length(x)-1); % xd的点数只有98个 
  > > plot(xd,df_c) 
  > > title('Derivative of fifth-deg. polynomial') 

以下的例子是针对数据组为离散型态，要注意的是原数据所代表的函数分布并无明显的 

折 
角，但是它的一阶微分经数值微分计算后的微分函数分布却有极大的曲折变化。 

  > > x=0:0.1:1; 
  > > y=[-.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.30 11.2]; 
  > > plot(x,y,'o',x,y) 
  > > title('y(x) data plot') 
  > > ylabel('y(x)'), xlabel('x') 
  > > dy=diff(y)./diff(x); 
  > > xd=x(1:length(x)-1); 
  > > plot(xd,dy) 
  > > title('Approximate derivative using diff') 
  > > ylabel('dy/dx'), xlabel('x') 

---------------------------------------------------------------------------- 

-- 
工程问题挑战：载具性能 
在此章我们简单的说明什么是常微分方程式，以及如何来应用数值方法求解方程式。 
在这章中提到的工程问题挑战「载具性能」，所谓载具(vehicle) 是泛指汽车、飞机、 

太 
空船。而这里谈的是飞机推进动力科技的一项进展，一种新式的无导管风扇(Unducted 
Fan, UDF) 的涡轮螺旋桨引擎(turboprop engine)。UDF的发展起源于1980年代中期，它 


利用新的材料、创新的叶片外型和更高叶片转速，可使螺旋桨客机飞行时速和喷射客机 

 ( 
以涡轮风扇为动力)一样快，而又更省油。 
---------------------------------------------------------------------------- 

-- 
-- 
第十章 常微分方程式 
10.1 微分方程式 
10.2 阮奇-库达方法 
10.3 范例问题：飞机发动机的加速性能分析 
10.4 高阶常微分方程式 
---------------------------------------------------------------------------- 

-- 
10.1 微分方程式 
---------------------------------------------------------------------------- 

-- 
-- 
一个一阶常微分方程式 (ordinary differential equation, ODE) 可以下式表示 
其中x 为独立变数而 y 是 x 的函数，我们就是要求解什么函数 y(x) 能满足上述 ODE 

。 
以下有几个一阶 ODE 的例子： 
除了上述的已知ODE外，还须有起始条件y0=y(x0)才能解方程式，即是在x=x0时，y(x)= 

y0 
。上述各个方程式的解析解 (analytical solution) 如下： 
以数值方法求解上述的 ODE 的问题，可以转换为在已知 y(a) 的函数值而要计算y(b)， 


依据泰勒序数对 y(b) 做展开 
其中 b=a+h。一个一阶的泰勒序数近似式为 
而二阶的泰勒序数近似式为 
MATLAB 所依据解ODE 的数值方法就是利用像上述的二阶及更高阶的三、四、五阶泰勒序 


数近似式来计算 f(b)。 
---------------------------------------------------------------------------- 

-- 
10.2 阮奇-库达方法 
---------------------------------------------------------------------------- 

-- 
-- 
阮奇-库达 (Runge-Kutta) 方法是最通用的解 ODE 的方法，它可以依计算精确度的要求 


有低阶到高阶的各个计算式，举例来说，一阶阮奇-库达法为 
其实上式即是一阶的泰勒序数近似式，只不过令  。因为一阶阮奇-库达法 的精确度太 

低 
，所以在真正解ODE 时，最少须用二阶以上的方法。 
MATLAB应用阮奇-库达法的函数有ode23, ode45，其中ode23是同时以二阶及三阶阮奇-库 


达法求解，而ode45 则是以四阶及五阶阮奇-库达法求解。其语法为ode23('dy',x0,xn, 

y0 
)，其中 dy 为ODE中的等式右边的函数（如之 前介绍的），x0, xn 是要解ODE的区间 
[x0, xn] 的二个端点，y0是起始值 (y0=y(x0))。而ode45的语法与ode23相同。 
先前的四个 ODE 的解法如下： 
例一、要在区间 [2, 4] 解以下的 ODE： 
% m-function, g1.m 
function dy=g1(x,y) 
dy=3*x.^2; 

  > > [x,num_y]=ode23('g1',2,4,0.5); 
  > > anl_y=x.^3-7.5; 
  > > plot(x,num_y,x,anl_y,'o') 
  > > title('Solution of g1') 
  > > xlabel('x'), ylabel('y=f(x)'), grid 

例二、要在区间 [0, 5] 解以下的 ODE： 
% m-function, g2.m 
function dy=g2(x,y) 
dy=-0.131*y; 

  > > [x,num_y]=ode23('g2',0,5,4); 
  > > anl_y=4*exp(-0.131*x); 
  > > plot(x,num_y,x,anl_y,'o') 
  > > title('Solution of g2') 
  > > xlabel('x'), ylabel('y=f(x)'), grid 

例三、要在区间 [0, 2] 解以下的 ODE： 
% m-function, g3.m 
function dy=g3(x,y) 
dy=2*x*cos(y)^2; 

  > > [x,num_y]=ode23('g3',0,2,pi/4); 
  > > anl_y=atan(x.*x+1); 
  > > plot(x,num_y,x,anl_y,'o') 
  > > title('Solution of g3') 
  > > xlabel('x'), ylabel('y=f(x)'), grid 

例四、要在区间 [0, 3] 解以下的 ODE： 
% m-function, g4.m 
function dy=g4(x,y) 
dy=3*y+exp(2*x); 

  > > [x,num_y]=ode23('g4',0,3,3); 
  > > anl_y=4*exp(3*x)-exp(2*x); 
  > > plot(x,num_y,x,anl_y,'o') 
  > > title('Solution of g4') 
  > > xlabel('x'), ylabel('y=f(x)'), grid 

如果将上述方法改以 ode45 计算，你可能无法察觉出其与ode23的解之间的差异，那是 

因 
为我们选的 ODE 代表的函数分布变化平缓，所以高阶方法就显现不出其优点。不过以二 


者在计算的误差上做比较，ode45 的误差量级会比 ode23要小。以下是几个例子： 
% m-function, g1.m 
function dy=g1(x,y) 
dy=3*x.^2; 
% m-file, odes1.m 
% Solve an ode using ode23 and ode45 
clg 
[x1,num_y1]=ode23('g1',2,4,0.5); 
anl_y1=x1.^3-7.5; 
error_1=abs(anl_y1-num_y1)./abs(anl_y1); % ode23 的计算误差 
[x2,num_y2]=ode45('g1',2,4,0.5); 
anl_y2=x2.^3-7.5; % 注意 x2 个数与 x1 不一定相同 
error_2=abs(anl_y2-num_y2)./abs(anl_y2); % ode45 的计算误差 
hold on 
subplot(2,2,1) 
plot(x1,num_y1,x1,anl_y1,'o') 
title('ODE23 solution'), ylabel('y') 
subplot(2,2,2) 
plot(x1,error_y1) % 注意二种方法的误差 
title('ODE23 error'), ylabel('y') % ode23 的误差的量级为 1.e-16 
subplot(2,2,3) 
plot(x2,num_y2,x2,anl_y2,'o') 
title('ODE45 solution'), ylabel('y') 
subplot(2,2,4) 
plot(x1,error_y2) 
title('ODE45 error'), ylabel('y') % ode45 的解没有误差 
hold off 
另一个例子： 
% m-function, g5.m 
function dy=g5(x,y) 
dy=-y+2*cos(x); 
% m-file, odes1.m 
% Solve an ode using ode23 and ode45 
clg 
[x1,num_y1]=ode23('g5',0,5,1); 
anl_y1=sin(x1)+cos(x1); 
error_1=abs(anl_y1-num_y1)./abs(anl_y1); 
[x2,num_y2]=ode45('g5',0,5,1); 
anl_y2=sin(x2)+cos(x2); 
error_2=abs(anl_y2-num_y2)./abs(anl_y2); 
hold on 
subplot(2,2,1) 
plot(x1,num_y1,x1,anl_y1,'o') 
title('ODE23 solution'), ylabel('y') 
subplot(2,2,2) 
plot(x1,error_y1) % 注意二种方法的误差 
title('ODE23 error'), ylabel('y') % ode23 的误差的量级为 1.e-4 
subplot(2,2,3) 
plot(x2,num_y2,x2,anl_y2,'o') 
title('ODE45 solution'), ylabel('y') 
subplot(2,2,4) 
plot(x1,error_y2) 
title('ODE45 error'), ylabel('y') % ode45 的误差的量级为 1.e-6 
hold off 
---------------------------------------------------------------------------- 

-- 
10.4 高阶常微分方程式 
---------------------------------------------------------------------------- 

-- 
-- 
一个高阶常微分方程式可以利用变数改变(change of variables) 方式改写成个一组联 

立 
的一阶常微分方程式。例如以下的 n 阶方程式 
我们先定义 n 个新的变数来取代上式中的 
将上述的新变数代入原 ODE，连同这些新变数即构成一组联立的一阶常微分方程式 
我们接者以一个二阶 ODE 为例说明上述的过程 
我们须要定义 
将上述的二变数代入原 ODE，即构成一组联立的一阶常微分方程式 
以下即是上述二阶 ODE 的解法： 
function u_prime =eqns2(x,u) 
u_prime(1) = u(1)*(1-u(2)^2) - u(2); 
u_prime(2) = u(1); 
initial = [0 0.25]; 
[x,num_y] = ode23('eqns2',0,20,initial); 
subplot(2,1,1), plot(x,num_y(:,1)) 
title('1st derivative of y'), xlabel('x'), grid 
subplot(2,1,2), plot(x,num_y(:,2)) 
title('y'), xlabel('x'), grid 
---------------------------------------------------------------------------- 

-- 
11.1.1 符号表示式 
---------------------------------------------------------------------------- 

-- 
-- 
在MATLAB中是将一符号表示式储存唯一字串 (character string)，即是以二个单引号之 


内的表示式来定义其为一符号式，例如 'tan(y/x)', 'x^3 - 2*x^2 + 3', 
'1/(cos(angle)+2)' 的三个式子。 
在一符号表示式中，需要定义所谓的独立变数。如果未曾事先指定何者为独立变数， 
MATLAB会自行决定。而它所决定变数的原则如下：它会挑选一个除了i和j之外而在字母 

上 
最接近x的小写字元；如果在式子中并无上述字元，则x会被视为预设的独立变数。函数 


symvar(s) 可以用已决定何者为独立变数。请看以下的例子： 
expression S symvar(S) 
'tan(y/x)' x 
'x^3-2*x^2+3' x 
'1/(cos(angle)+2)' x 
'3*a*b-6' b 
MATLAB提供了一个函数ezplot 可以画单变数的符号式，其预设的独立变数的范围是。它 


的语法为 ezplot(S), S代表符号变数；另一各相关语法 ezplot(S,[xmin,xmax])，则是 


设定独立变数的范围xmin到xmax。