matlab第八讲数值计算.ppt

Matlab第八讲数值计算 赵瑞 数值计算 线性方程组的解函数的零点函数的极值点数值积分曲线的多项式拟合 一 线性方程组的解 线性方程组AX b A为可逆矩阵 求X的命令 X inv A b b输入成列向量 X返回列向量不推荐使用逆矩阵法2 推荐使用 除法 精度比逆矩阵法高X A b A左除b 相当于A 1 b b输入成列向量 X返回列向量3 X b A b右除A 相当于b A 1 b输入成行向量 X返回行向量例 A 1 1 2 1 3 1 2 1 3 1 4 1 3 1 4 1 5 b 1 2 3 X1 inv A bX2 A bb b X3 b A 二 函数的零点 对于任意函数f x 0来说 它可能有零点 也可能没有零点 可能只有一个零点 也可能有多个甚至无数个零点 因此 很难说出一个通用解法 一般说来 零点的数值计算过程是 先猜测一个初始零点或该零点所在的区间 然后通过一些计算 使猜测值精确化 利用fzero指令求一元函数零点 x y fzero fun x0 fzero是根据函数是否穿越横轴来决定零点的 因此本程序无法确定函数曲线仅触及横轴和不穿越的零点 例如 y x2的零点 不穿越横轴 fun表示函数 数据类型是字符串x0是表示零点初始猜测 x0可以是标量或二元向量 当x0取标量时 该指令将在它两侧寻找一个与之最靠近的零点 当x0取向量 a b 时 该指令将在区间 a b 内寻找一个零点 返回值x表示零点 y表示在这点的函数值注意 这个命令只能求x0附近的一个零点 例 求函数f x sin2x e 0 1x 0 5 x 的零点 1 首先用字符串表示函数 matlab命令 f sin x 2 exp 0 1 x 0 5 abs x 2 作图法观察函数零点分布 matlab命令 x 10 0 01 10 y sin x 2 exp 0 1 x 0 5 abs x plot x y holdon 接着画 下面画出x轴xx 10 0 01 10 yy zeros size xx plot xx yy r 观察f x 与x轴的交点 大致可以看出有四个穿越x轴的交点 具体位置不是很清楚 所以希望图像能够局部放大 以便看得更加清楚 3 利用zoom和ginput指令获得零点的初始近似值 matlab命令 zoomon 在MATLAB指令窗中运行 获局部放大图 或者直接在图形窗口单击放大镜图标即可 局部放大之后的图形 xx yy ginput 5 在MATLAB指令窗中运行 用鼠标获5个零点猜测值 zoomoffholdoff 4 求靠近xx各点的精确零点 matlab命令 x1 y1 fzero f xx 1 求在xx 1 附近的精确零点 x2 y2 fzero f xx 2 求在xx 2 附近的精确零点 x3 y3 fzero f xx 3 求在xx 3 附近的精确零点 x4 y4 fzero f xx 4 求在xx 4 附近的精确零点 x5 y5 fzero f xx 5 求在xx 5 附近的精确零点 运行结果x1 2 0074y1 2 2204e 016x2 0 5198y2 0 x3 0 5198y3 0 x4 0 5993y4 5 5511e 017x5 1 6738y5 1 1102e 016 5 求在xx 3 附近的精确零点从理论分析可知 函数f x sin2x e 0 1x 0 5 x 0是函数的一个零点 但即便是以十分靠近该零点的值为搜索的初始值 也找不到 而却找到了另一个零点 原因是曲线没有穿越横轴 三 求函数极值点 Matlab函数中只有处理极小值的指令 求的是局部极值1 一元函数极值 x fval fminbnd fun x1 x2 x1 x2表示被研究区间的左右边界 x返回极小值点 fval返回该点的函数值 2 多元函数极值 x fval fminsearch fun x0 单纯形法求多元函数极值点指令 x0为初始值 x fval fminunc fun x0 拟牛顿法求多元函数极值点指令 x0为初始值 fun均是用字符串表示函数 x fval fminbnd sin x pi pi 求sinx在 的极小值点x 1 5708fval 1 0000 x fval fminsearch x 1 1 2 x 2 x 1 2 2 0 0 x 1 00001 0000fval 2 6131e 010 x fval fminunc x 1 2 x 2 2 x 1 x 2 2 x 1 4 x 2 2 2 x 02fval 4 四 数值积分 梯形方法trapz辛普森公式quad科特斯公式quad8精度依次提高 梯形方法matlab命令 x eps 0 001 1 eps表示什么 y sin x x 注意 trapz x y ans 0 9461 由于辛普森公式quad 科特斯公式quad8这二种方法依据的积分法不同于梯形法 因此它们的语法就和trapz不同 其语法为quad function a b quad8语法相同 其中function是一已定义函数的名称 如sin cos sqrt log等或自己定义的函数 而a b是积分的下限和上限 和trapz比较 quad quad8不同之处在于这二者类似解析式的积分式 只须设定上下限及定义要积分的函数 而trapz则是针对离散点型态的数据做积分 quad sin x x 0 1 Warning Dividebyzero Ininlineevalat13Ininline subsrefat25Inquadat62ans 0 9461 quad sin x x eps 1 ans 0 9461 我们看一简单积分式以下为MATLAB的程序formatlongkq quad sqrt 0 1 kq 0 66665956927202kq8 quad8 sqrt 0 1 kq8 0 66666534581698 五 曲线的多项式拟合 1 多项式拟合已知一些离散的数据点x x0 x1 xm 和对应的函数值y y0 y1 ym 用n次多项式拟合这组数据 p s polyfit x y n Fitpolynomialtodata 拟合的n次多项式为p x anxn an 1xn 1 a1x a0p返回的是一个向量p an an 1 a1 a0 s是一个构架数组 其中s normr表示误差 例 已知一组实验数据如下 求它的拟合曲线首先画图观察x 12345 y 44 5688 5 plot x y r 点的形状比较接近一条直线 所以用一次多项式拟合 p s polyfit x y 1 p 1 25002 4500s R 2x2double df 3normr 0 8216 下面画拟合曲线的图 首先取点xn 1 0 001 5 yn polyval p xn Evaluatepolynomial 若p an an 1 a1 a0 则此命令求多项式anxn an 1xn 1 a1x a0在xn各分量处的函数值plot xn yn b x y r 2 其他形式 1 已知一些离散的数据点x x0 x1 xm 和对应的函数值y y0 y1 ym 拟合成型如y aebx做变换lny bx lna 这是个一次多项式 可以用xi lnyi拟合 i 0 1 m 拟合的曲线为 yy bx lna用来拟合的数据 x x0 x1 xm t lny0 lny1 lnym 命令 x x0 x1 xm y y0 y1 ym yy log y p s polyfit x yy 1 此时p返回的是 blna 思考 如果想画y aebx的图像该怎么做 此时p返回的是 blna 思考 如果想画y aebx的图像该怎么做 答案 xn x0 步长 xm 取一些点 由于拟合的曲线是y aebx 计算在xn各分量处的函数值yn exp p 2 exp p 1 xn plot x y r xn yn b 画图legend 数据点 拟合曲线 最好加上图例标注 例 matlab命令 x 1 00 1 25 1 50 1 75 2 00 y 5 10 5 79 6 53 7 45 8 46 plot x y r 通过观察 确定拟合曲线形如y a exp b x 离散数据点的图像 通过观察确定拟合曲线形如y a exp b x yy log y p s polyfit x yy 1 p 1 b p 2 lna 下面画图y a exp b x xn 1 0 001 2 yn exp p 2 exp p 1 xn holdonplot xn yn legend 数据点 拟合曲线 holdoff 运行结果 p 0 50571 1225s R 2x2double df 3normr 0 0052 作业 衡量一段时间 t1 t2 的用药疗效 通常看这段时间cd4细胞的净增率 t2时刻的cd4细胞浓度 t1时刻的cd4细胞浓度 t1时刻的cd4细胞浓度对附件2中的数据进行分析 要求 1 只考虑测量了6次的病人 并且各次cd4细胞浓度均不为零2 对每种疗法按年龄段分类 两个年龄段 40岁以下 包括40 40岁以上 不包括40 3 观察疗效并且预测以后的疗效例