第六章 计算方法  非线性方程求解 多项式插值与曲线拟合 数值微分与数值积分 求常微分方程数值解命令

已知函数 f(x),求 x 使 f(x)=0 —— 非线性方程 例 x = tan x x y (k = 1,2,·······) f=inline('x-tan(x)') x=fzero(f,3.5) x = 4.4934 求一元函数零点 fzero(fun, x0) X = fzero(fun, x0)寻找fun在x0附近的零点

r x 例 x 3 – 30 x 2 + 2552 =0 p=[1 -30 0 2552]; roots(p) 半径r =10 cm的球体,密度 =0.638.浸入水深度 x=? ans = 26.3146 11.8615 -8.1761 15/18

求多项式零点: roots(p); 例2 求 x3 –8x2 +6x– 30=0的解 P=[1 -8 6 -30]; r=roots(P) r = 7.7260 0.1370 + 1.9658i 0.1370 - 1.9658i 例1.求方程 x3 – 4x + 5 = 0 的解 P=[1 0 -4 5]; roots(P) ans = -2.4567 1.2283 + 0.7256i 1.2283 - 0.7256i

P10(t) f(t) f(x) 多项式插值原理 已知 f(x) 在点 xi 上的函数值 yi=f(xi), (i=0,1,2,···,n) 求多项式 P(x)=a0 + a1x +···+ anxn 满足: P(xk)= yk (k = 0,1,…,n) P10(t) f(t) f(x)

zi=interp2(x, y, z, xi, yi, ‘method’) 一维插值: yi = interp1(x, y, xi, ‘method ’) x = 0:10; y = sin(x); xi = 0:.25:10; yi = interp1(x,y,xi); plot(x,y,'o',xi,yi) method nearest 最近点插值 linear 线性插值 spline 样条插值 cubic 立方插值 二维插值 zi=interp2(x, y, z, xi, yi, ‘method’) 三维插值 vi = interp3(x,y,z,v,xi,yi,zi, ‘method’) 17/18

例1.正弦曲线数据插值试验

离散数据的多项式拟合 x x1 x2 ·········· xm f(x) y1 y2 ·········· ym 使得 已知数据表 求拟合函数:  (x) = a0 + a1x + ······+ anxn 使得 达到最小

P(x)=a0xn + a1xn-1 + ······ + an-1x + an 多项式拟合命令 P=polyfit(X,Y,n) 求出(最小二乘意义下)n次拟合多项式 P(x)=a0xn + a1xn-1 + ······ + an-1x + an 计算结果为系数: P=[ a0, a1, ······ , an-1, an ] 多项式求值命令: y1=polyval(P, x) 其中,P是n次多项式的系数,x是自变量的值,y1是多项式在x处的值 17/18

定积分与积分和式 Sn 5.2908 5.1044 4.9835 ······ 4.8999 右矩形和 h 1 0.5 0.2 ······ 4/16

梯形数值积分公式 左矩形 梯形 右矩形 4.4429 4.8669 5.2908 4.6804 4.8924 5.1044 4.8139 4.8987 4.9835 4.8572 4.8996 4.9420 5/16

梯形数值积分命令 trapz() clear x=0:pi/100:pi; y=sin(x); trapz(x,y) ans = 1.9998 clear x=sort(rand(1,101)*pi); y=sin(x); trapz(x,y) ans = 1.9981 rand(1,101)产生101个均匀随机数,每个数都介于0-1之间

辛卜生求数值积分命令quad() clear fun=inline('1./(x.^3-2*x-5)') ezplot(fun,[0,2]) [q,n]=quad(fun,0,2) q = -0.4605 n = 53

二重积分dblquad()与三重积分 fun=inline('y*sin(x)+x*cos(y)') Q=dblquad(fun,pi,2*pi,0,pi) Q = -9.8698 [x,y]=meshgrid(pi:.1:2*pi,0:.1:pi); z=fun(x,y); mesh(x,y,z)

物理量在空间的分布情况其及随时间变化的规律。 研究具有实际应用背景的数学函数 物理量在空间的分布情况其及随时间变化的规律。 物理现象  物理定律  微分方程  数值解 建模型: 物理学家、化学家、生物学家、社会学家 2/18

 y(t)=y0 e r t 马尔萨斯(Malthus)人口模型 某时刻人口数为y0, 在以后某一时间 t 人口数y 的变化率正比于y Thomas Malthus (1766-1834) MATLAB符号求解方法 y=dsolve('Dy=r*y','t') y =C1*exp(r*t) 弹性力学问题 单摆问题 天体力学二体问题 电磁场计算············ 令 r = a, y0 = exp( b )  y(t) = exp( at + b) a = ?, b = ?  人口数据拟合实验 3/18

MATLAB求常微分方程数值解 简例 f=inline('y-x.*y.^2'); [x,y]=ode23(f,[0,2],1) 0 1.0000 0.0800 1.0796 0.2621 1.2485 0.4437 1.3755 0.6437 1.4402 0.8437 1.4209 1.0437 1.3374 1.2437 1.2196 1.4437 1.0923 1.6437 0.9709 1.8437 0.8620 2.0000 0.7870 f=inline('y-x.*y.^2'); [x,y]=ode23(f,[0,2],1) 4/18

ode23(),ode45() 求常微分方程初值问题数值解函数 解一阶常微分方程组初值问题命令格式 [T,y] = ode23(' F ',Tspan,y0) 其中, F是ODE文件, 表示 y' = F(t,y)中右端函数 Tspan = [t0 Tfinal] —— 求解区域; y0 —— 初始条件 注: 函数F(t,y) 必须返回列向量. 数值解 y 的每一行对应于列向量T中的每一时刻 5/18

捕食者和被捕食者模型 ( 0≤ t ≤ 15 ) y1(0)=20, y2(0)=20 数会急增，数学家将鱼分成捕食者和被捕食者，分别 生物学家研究发现，战争期间捕鱼量减少,食肉鱼总 数会急增，数学家将鱼分成捕食者和被捕食者，分别 以y2(t) , y1(t)表示它们在时刻t 的总数。假定被捕食者 (食用鱼)有充分的食物。 ( 0≤ t ≤ 15 ) y1(0)=20, y2(0)=20 6/18

function yp = lotka(t,y) % predator-prey model. MATLAB求解: (1)定义函数; (2)调用ode23 函数文件,文件名:lotka.m function yp = lotka(t,y) % predator-prey model. yp=diag([1 - .01*y(2), -1+.02*y(1)])*y; t0=0;tf=15; y0=[20 20]; [t,y]=ode23(‘lotka’,[t0,tf],y0); plot(t,y(:,1),t,y(:,2),'r') —— 捕食者y2 —— 被捕食者y1 7/18

洛伦兹模型与“蝴蝶效应” x(0)=0，y(0)=0，z(0)=0.01 function z=flo(t,y) Y-X平面 x(0)=0，y(0)=0，z(0)=0.01  = 8/3  = 10  =28 function z=flo(t,y) A=[-8/3 0 y(2); 0 -10 10; -y(2) 28 -1]; z=A*y; y0=[0;0;0.01]; [x,y]=ode23(‘flo’,[0, 80],y0); plot(y(:,2),y(:,1)) 8/18