MATLAB 语言程序设计基础(2)

%第六章微分方程问题的解法% 微分方程的解析解方法% 常微分方程问题的数值解法% 微分方程问题算法概述% 四阶定步长Runge-Kutta算法及MATLAB 实现% 一阶微分方程组的数值解% 微分方程转换% 特殊微分方程的数值解% 边值问题的计算机求解% 偏微分方程的解% 6.1 微分方程的解析解方法% y=dsolve(f1, f2, …, fm ,'x')% syms t; u=exp(-5*t)*cos(2*t+1)+5;% uu=5*diff(u,t,2)+4*diff(u,t)+2*u% syms t y;% y=dsolve(['D4y+10*D3y+35*D2y+50*Dy+24*y='...%'87*exp(-5*t)*cos(2*t+1)+92*exp(-5*t)*sin(2*t+1)+10'],'y(0)=3','Dy(0)=2','D2y(0)=0','D3y(0)=0')% [x,y]=dsolve('D2x+2*Dx=x+2*y-exp(-t)',...% 'Dy=4*x+3*y+4*exp(-t)')% syms t x;% x=dsolve('Dx=x*(1-x^2)+1')% Warning: Explicit solution could not be found; implicit solution returned.% > In D:\MATLAB6p5\toolbox\symbolic\dsolve.m at line 292% x =% t-Int(1/(a-a^3+1),a=``..x)+C1=0% 故只有部分非线性微分方程有解析解。

% 6.2 微分方程问题的数值解法% 6.2.1 微方程问题算法概述%Euler算法%%function [outx,outy]=MyEuler(fun,x0,xt,y0,PointNum)% fun 表示f(x,y); x0,xt:自变量的初值和终值;% y0:函数在x0处的值,其可以为向量形式;% PointNum表示自变量在[x0,xt]上取的点数% if nargin<5 |PointNum<=0 %PointNum 默认值为100% PointNum=100;% end% if nargin<4 % y0默认值为0% end% h=(xt-x0)/PointNum; %计算步长h% x=x0+[0:PointNum]'*h; %自变量数组% y(1,:) = y0(:)'; %将输入存为行向量，输入为列向量形式% for k = 1:PointNum% f=feval(fun,x(k),y(k,:)); f=f(:)'; %计算f(x,y)在每个迭代点的值% y(k + 1,:) =y(k,:) +h*f; %对于所取的点x迭代计算y值% end% outy=y; outx=x;% plot(x,y) %画出方程解的函数图% end% [x1,y1]=MyEuler('myfun',0,2*pi,1,16);% myfun=incline'Dy=sin(x)+y'% function [Xout,Yout]=MyEulerPro(fun,x0,xt,y0,PointNumber)% %MyEulerPro 用改进的欧拉法解微分方程% if nargin<5 | PointNumber<=0% %PointNumer默认值为100% PointNumer=100;% end% if nargin<4 %y0默认值为0% y0=0;% end% h=(xt-x0)/PointNumber;% %计算所取的两离% 散点之间的距离% x=x0+[0:PointNumber]'*h;% %表示出离散的自变量x% y(1,:)=y0(:)';% for i=1:PointNumber %迭代计算过程% f1=h*feval(fun,x(i),y(i,:)); f1=f1(:)';% f2=h*feval(fun,x(i+1),y(i,:)+f1); f2=f2(:)';% y(i+1,:)=y(i,:)+1/2*(f1+f2);% end% Xout=x;Yout=y;%6.2.2 四阶定步长Runge-Kutta算法及MATLAB 实现%6.2.3 一阶微分方程组的数值解%6.2.3.1 四阶五级Runge-Kutta-Felhberg算法%6.2.3.2 基于MATLAB 的微分方程% 求解函数% 格式1：直接求解% [t,x]=ode45(Fun,[t0,tf],x0)% 格式2：带有控制参数% [t,x]=ode45(Fun,[t0,tf],x0，options)% 格式3：带有附加参数% *t,x+=ode45(Fun,*t0,tf+,x0,options,p1,p2,…)% [t0,tf]求解区间，x0初值问题的初始状态变量。

% 描述需要求解的微分方程组：% 不需附加变量的格式% function xd=funname(t,x)% 可以使用附加变量% function xd=funname(t,x,flag,p1,p2,…)% t是时间变量或自变量（必须给），x为状态向量，xd为返回状态向量的导数. % flag用来控制求解过程,指定初值，即使初值不用指定，也必须有该变量占位。

% 修改变量：options唯一结构体变量，用odeset( )修改。

% options=odeset(‘RelTol’,1e-7);% options= odeset; options. RelTol= 1e-7;%6.2.3.3 MATLAB 下带有附加参数的微分方程求解% 6.2.4 微分方程转换% 6.2.4.1 单个高阶常微分方程处理方法% 6.2.4.2 高阶常微分方程组的变换方法%% 6.3特殊微分方程的数值解% 6.3.1 刚性微分方程的求解% MATLAB采用求解函数ode15s()，该函数的调用格式和ode45()完全一致。

% *t,x+=ode15s(Fun,*t0,tf+,x0,options,p1,p2,…)%% 6.3.2 隐式微分方程求解% 6.3.3 微分代数方程求解% 6.3.4延迟微分方程求解%% 6.4边值问题的计算机求解% 6.4.1 边值问题的打靶算法% 6.4.2 线性微分方程的有限差分算法%% 6.5 偏微分方程求解入门% 6.5.1 偏微分方程组求解% pdepe()函数可求解% 函数描述% [c,f,s]=pdefun(x,t,u,ux)% 边界条件函数描述% [pa,qa,pb,qb]=pdebc(x,t,u,ux)% 初值条件函数描述% u0=pdeic(x)% 求解偏微分方程% sol=pdepe(m,@pdefun,@pdeic,@pdebc,x,t)%% 6.5.2 二阶偏微分方程的数学描述% 椭圆型；抛物线型；双曲线型；特征值型偏微分方程%% 6.5.3 偏微分方程的求解界面应用简介% 6.5.3.1 偏微分方程求解程序概述% 启动偏微分方程求解界面% 在MATLAB 下键入pdetool% 该界面分为四个部分：菜单系统；工具栏；集合编辑；求解区域% 6.5.3.2 偏微分方程求解区域绘制% 1）用工具栏中的椭圆、矩形等绘制一些区域。

% 2）在集合编辑栏中修改其内容。

如（R1＋E1＋E2）－E3% 3）单击工具栏中按纽可得求解边界。

% 4）选择Boundary-Remove All Subdomain Borders菜单项，消除相邻区域中间的分隔线。

% 5）单击按纽可将求解区域用三角形划分成网格。

可用按纽加密。

%% 6.5.3.3 偏微分方程边界条件描述% 选择Boundary-Specify Boundary Conditions菜单% 狄利克雷条件，诺伊曼条件。

% 6.5.3.4 偏微分方程求解举例% 6.5.3.5 函数参数的偏微分方程求解% 第七章代数方程与最优化问题的求解% 代数方程的求解% 无约束最优化问题的计算机求解% 有约束最优化问题的计算机求解% 整数规划问题的计算机求解% 7.1代数方程的求解% 7.1.1 代数方程的图解法% 7.1.2 多项式型方程的准解析解法% 由于方程阶次可能太高，不存在解析解。

然而，% 可利用MATLAB的符号工具箱得出原始问题的高精度数值解，故称之为准解析解。

% 一般多项式方程的根可为实数，可为复数。

MATLAB符号工具箱中的solve( )函数。

% 最简调用格式：S=solve(eqn1,eqn2,…,eqnn)% (返回一个结构题型变量S，如S.x表示方程的根。

）% 直接得出根：（变量返回到MATLAB工作空间）% *x,…+=solve(eqn1,eqn2,…,eqnn)% 同上，并指定变量% [x,…]=solve(eqn1,eqn2,…,eqnn，’x,…')% 7.1.3 一般非线性方程数值解% 非线性方程的标准形式为f(x)=0% 函数fzero% 格式x = fzero (fun,x0) %用fun定义表达式f(x)，x0为初始解。

% x = fzero (fun,x0,options)% *x,fval+=fzero(…) %fval=f(x)% *x,fval,exitflag+=fzero(…)% *x,fval,exitflag,output+=fzero(…)% 说明该函数采用数值解求方程f(x)=0的根。

%% 7.1.4 一般非线性方程组数值解% 格式：最简求解语句% x=fsolve(Fun, x0)% 一般求解语句% *x, f, flag, out+=fsolve(Fun, x0,opt, p1, p2,…)% 若返回的flag 大于0，则表示求解成功，否则求解出现问题，% opt 求解控制参数，结构体数据。

% 获得默认的常用变量opt=optimset;% 用这两种方法修改参数（解误差控制量）% opt.Tolx=1e-10;或set(opt.‘Tolx’, 1e-10)% 可先用用图解法选取初值，再调用fsolve( )函数数值计算% 7.2无约束最优化问题求解% 7.2.1 解析解法和图解法%% 单变量函数求最小值% 函数fminbnd (最值可能在端点，需要考虑）% 格式x = fminbnd(fun,x1,x2)% x = fminbnd(fun,x1,x2,options)% [x,fval] = fminbnd(…) % fval为目标函数的最小值% *x,fval,exitflag+ = fminbnd(…)% exitflag为终止迭代的条件,若exitflag>0，收敛于x，exitflag=0，% 表示超过函数估计值或迭代的最大数字，exitflag<0表示函数不收敛于x；% *x,fval,exitflag,output+ = fminbnd(…)% output为优化信息,若output=iterations表示迭代次数，% output=funccount表示函数赋值次数，output=algorithm表示所使用的算法%% 7.2.2 基于MATLAB的数值解法%求解无约束最优化的函数：fminsearch()% x=fminunc(Fun,x0) %最简求解语句% [x,f,flag,out]=fminunc(Fun,x0,opt,p1,p2,...)% 一般求解格式% 比较可知fminunc()函数效率高于fminsearch()函数，但% 当所选函数高度不连续时，使用fminsearch效果较好。