华科计算方法上机实验报告.doc

1 计算方法计算方法 实验报告实验报告 专业及班级专业及班级 姓名 姓名 学号 学号 日期 日期 20132013年年6 6月月1212日日 一 方程求根一 方程求根 1 用牛顿迭代法求解下列方程的正根 1 log 1 x x 2 0 x 7 468818e 001 求解代码 牛顿迭代法的函数 M 文件 function x Newt n f name x0 x x0 xb x 1 k 0 n 50 del x 0 01 while abs x xb 0 000001 k k 1 xb x if k n break end y feval f name x y driv feval f name x del x y del x x xb y y driv fprintf k 3 0f x 12 5e y 12 5e yd 12 5e n k x y y driv end fprintf n Final answer 12 6e n x 待求方程的函数 M 文件 function y eqn 1 x y log x 1 x 2 运行结果 Newt n eqn 1 1 k 1 x 7 96954e 001 y 3 06853e 001 yd 1 51125e 000 k 2 x 7 50200e 001 y 4 90424e 002 yd 1 04895e 000 k 3 x 7 46936e 001 y 3 07003e 003 yd 9 40663e 001 k 4 x 7 46882e 001 y 5 03365e 005 yd 9 33074e 001 k 5 x 7 46882e 001 y 6 30906e 007 yd 9 32949e 001 Final answer 7 468818e 001 2 e x 5x 0 x 6 052671e 001 待解方程 function y eqn 1 x 2 y exp x 5 x 2 运行结果 Newt n eqn 2 2 k 1 x 1 00102e 000 y 1 26109e 001 yd 1 26239e 001 k 2 x 6 88531e 001 y 2 28918e 000 yd 7 32552e 000 k 3 x 6 11607e 001 y 3 79584e 001 yd 4 93453e 000 k 4 x 6 05365e 001 y 2 69226e 002 yd 4 31343e 000 k 5 x 6 05268e 001 y 4 13312e 004 yd 4 26254e 000 k 6 x 6 05267e 001 y 3 99547e 006 yd 4 26175e 000 Final answer 6 052671e 001 ans 0 6053 3 x 3 2x 1 0 x 4 533977e 001 待解方程 function y eqn 1 x y x 3 2 x 1 运行结果 Newt n eqn 3 3 k 1 x 1 89997e 000 y 3 20000e 001 yd 2 90901e 001 k 2 x 1 15047e 000 y 9 65861e 000 yd 1 28868e 001 k 3 x 6 80277e 001 y 2 82368e 000 yd 6 00536e 000 k 4 x 4 82154e 001 y 6 75369e 001 yd 3 40884e 000 k 5 x 4 53984e 001 y 7 63945e 002 yd 2 71198e 000 k 6 x 4 53401e 001 y 1 53570e 003 yd 2 63202e 000 k 7 x 4 53398e 001 y 8 46834e 006 yd 2 63042e 000 k 8 x 4 53398e 001 y 4 41262e 008 yd 2 63041e 000 Final answer 4 533977e 001 ans 0 4534 4 x 2 0 5 x 0 x 2 0000 待解方程 function y eqn 1 x y x 2 0 5 x 运行结果 Newt n eqn 4 1 k 1 x 2 02879e 000 y 7 32051e 001 yd 7 11565e 001 k 2 x 2 00002e 000 y 2 16052e 002 yd 7 51049e 001 k 3 x 2 00000e 000 y 1 72788e 005 yd 7 50157e 001 k 4 x 2 00000e 000 y 3 60276e 009 yd 7 50156e 001 Final answer 2 000000e 000 ans 2 0000 3 2 先用图解法确定初始点 然后再求方程 x sin x 1 0 x 0 10 的所有根 画图程序 clf hold on grid on x 0 0 01 10 y1 x sin x 1 y2 0 plot x y1 r plot x y2 g 运行结果 从图中可以得出该方程的几个根分别接近于 1 3 6 5 和 9 则可以进行下一步运算 待解方程 function y eqn 1 x y x sin x 1 运行结果 Newt n eqn 5 1 k 1 x 1 11463e 000 y 1 58529e 001 yd 1 38292e 000 k 2 x 1 11416e 000 y 6 61864e 004 yd 1 38810e 000 k 3 x 1 11416e 000 y 3 24115e 007 yd 1 38817e 000 Final answer 1 114157e 000 ans 1 1142 4 Newt n eqn 5 3 k 1 x 2 79702e 000 y 5 76640e 001 yd 2 84083e 000 k 2 x 2 77312e 000 y 5 51756e 002 yd 2 30892e 000 k 3 x 2 77261e 000 y 1 14804e 003 yd 2 24109e 000 k 4 x 2 77260e 000 y 7 70019e 006 yd 2 23963e 000 k 5 x 2 77260e 000 y 4 91882e 008 yd 2 23962e 000 Final answer 2 772605e 000 ans 2 7726 Newt n eqn 5 6 5 k 1 x 6 43934e 000 y 3 98280e 001 yd 6 56560e 000 k 2 x 6 43912e 000 y 1 44090e 003 yd 6 52127e 000 k 3 x 6 43912e 000 y 1 07492e 006 yd 6 52106e 000 Final answer 6 439117e 000 ans 6 4391 Newt n eqn 5 9 k 1 x 9 34662e 000 y 2 70907e 000 yd 7 81559e 000 k 2 x 9 31742e 000 y 2 70260e 001 yd 9 25348e 000 k 3 x 9 31724e 000 y 1 59253e 003 yd 9 17140e 000 k 4 x 9 31724e 000 y 2 61235e 006 yd 9 17089e 000 Final answer 9 317243e 000 ans 9 3172 由此可以解得方程在 0 10 之间的四个根分别为 1 1142 2 7726 6 4391 和 9 3172 二 线性方程组二 线性方程组 1 计算下列矩阵的逆矩阵 并验证之 1 运行过程 a 1 4 5 2 1 2 8 1 1 c inv a c 0 0400 0 0400 0 1200 0 5600 1 5600 0 3200 0 2400 1 2400 0 2800 d c a d 1 0000 0 0000 0 0000 0 0000 1 0000 0 0000 0 0000 0 1 0000 e a c e 1 0000 0 0000 0 0000 0 0000 1 0000 0 0000 0 0000 0 0000 1 0000 5 2 运行过程 a 0 5 1 1 6 3 8 9 5 c inv a c 0 5377 0 3208 0 0849 0 2736 0 0755 0 0094 0 3679 0 3774 0 0472 d c a d 1 0000 0 0000 0 0000 0 1 0000 0 0000 0 0 0000 1 0000 e a c e 1 0000 0 0000 0 0000 0 0000 1 0000 0 0000 0 0000 0 0000 1 0000 2 用两种方法求解下列线性方程组 调用 x A b 命令 1 x 2 5000 4 0000 3 5000 T 2 x 1 0000 2 6667 8 6667 T 运行过程 a 2 1 0 1 2 1 0 1 2 b 1 2 3 x a b x 2 5000 4 0000 3 5000 c 2 1 1 3 4 1 0 1 1 d 4 5 6 y c d y 1 0000 2 6667 8 6667 利用 LU 分解 1 运行过程 a 2 1 0 1 2 1 0 1 2 b 1 2 3 l u p lu a l 1 0000 0 0 0 5000 1 0000 0 0 0 6667 1 0000 u 6 2 0000 1 0000 0 0 1 5000 1 0000 0 0 1 3333 p 1 0 0 0 1 0 0 0 1 调用命令 y l p b 求得 y 1 0000 2 5000 4 6667 T 调用命令 x u y 求得 x 2 5000 4 0000 3 5000 T 2 运行过程 c 2 1 1 3 4 1 0 1 1 d 4 5 6 l u p lu c l 1 0000 0 0 0 6667 1 0000 0 0 0 6000 1 0000 u 3 0000 4 0000 1 0000 0 1 6667 0 3333 0 0 1 2000 p 0 1 0 1 0 0 0 0 1 调用命令 y l p b 求得 y 5 0000 7 3333 10 4000 T 调用命令 x u y 求得 x 1 0000 2 6667 8 6667 T 三 插值与拟合三 插值与拟合 1 lagrange 方程 M 文件 function fi lagrange x y xi fi zeros size xi npl length y for i 1 npl z ones size xi for j 1 npl if i j z z xi x j x i x j end end fi fi z y i end 运行过程与结果 x 0 0 4 0 8 1 2 y 1 0 1 491 2 225 3 320 xi 0 2 0 6 1 0 yi lagrange x y xi 7 yi 1 2225 1 8203 2 7200 所以得到三个插值点的 yi 分变为 1 2225 1 8203 和 2 7200 误差计算 k exp xi yi k 0 0011 0 0019 0 0017 所以误差分别为 0 0011 0 0019 和 0 0017 2 1 飞机降落曲线 y 1 1574074e 009 x 3 2 0833333e 005 x 2 运行过程 由题意可以得到 a1 a0 都为 0 则有 a 12000 3 12000 2 3 12000 12000 2 12000 b 1000 0 x a b x 1 0e 004 0 0000 0 2083 由于计算结果太小无法完全显示 故先保存为 txt 文档后再得到 a3 a2 的值 save Adat txt ascii 最终可以得到 a3 1 1574074e 009 a2 2 0833333e 005 a1 a0 0 于是方程解得为 y 1 1574074e 009 x 3 2 0833333e 005 x 2 2 绘制出飞机降落曲线 运行过程 clf x 0 0 1 20000 y 1 1574074e 009 x 3 2 0833333e 005 x 2 plot x y grid on 得到图像如下 8 3 用三次多项式拟合下面数据 并做出图形 运行代码 x 0 0 2 0 4 0 6 0 8 1 y 0 7 78 10 68 8 37 3 97 0 c polyfit x y 3 c 41 5625 101 6071 60 0768 0 1179 于是可以得到拟合方程 y 41 5625 x 3 101 6071 x 2 60 0768 x 0 1179 画图 clf x 2 0 001 2 y 41 5625 x 3 101 6071 x 2 60 0768 x 0 1179 plot x y grid on 得到图形 四 数值积分四 数值积分 1 用复合梯形求积法计算下列积分 取 n 2 4 8 16 1 梯形公式的函数 M 文件 function t trapezia a b n f name h b a n fa feval f name a fb feval f name b t1 0 for k 1 n 1 x k a k h t1 t1 feval f name x k end t h fa 2 t1 fb 2 9 原方程代码 function y func ts x y exp x 运行结果 t trapezia 0 1 2 func ts t 1 7539 t trapezia 0 1 4 func ts t 1 7272 t trapezia 0 1 8 func ts t 1 7205 t trapezia 0 1 16 func ts t 1 7188 所以可以得到 T2 1 7539 T4 1 7272 T8 1 7205 T16 1 7188 2 方程代码 function y func ts x y 1 2 x 运行结果 t trapezia 0 1 2 func ts t 0 4083 t trapezia 0 1 4 func ts t 0 4062 t trapezia 0 1 8 func ts t 0 4056 t trapezia 0 1 16 func ts t 0 4055 所以可以得到 T2 0 4083 T4 0 4062 T8 0 4056 T16 0 4055 2 用复合 Simpson 求积法计算下列积分 取 n 4 8 16 32 1 复合 Simpson 公式的函数 M 文件 function s simpson a b n f name h b a n for k 1 2 n x k a k h 2 end fa feval f name a fb feval f name b s1 0 s2 0 for k 1 n 10 s1 s1 feval f name x 2 k 1 s2 s2 feval f name x 2 k end s h 0 5 fa fb 2 s1 s2 3 待解方程 function y func ts x y 1 2 cos x 运算结果 s simpson 0 pi 4 func ts s 1 8138 s simpson 0 pi 8 func ts s 1 8138 s simpson 0 pi 16 func ts s 1 8138 s simpson 0 pi 32 func ts s 1 8138 所以可以得到 S4 1 8138 S8 1 8138 S16 1 8138 S32 1 8138 2 待解方程 function y func ts x y log 1 x x 求解过程 由于待解方程在 x 0 处没有意义 所以将积分起点设定为很接近 0 的一个正数 s simpson 0 0000000000000000001 2 4 func ts s 1 3534 s simpson 0 0000000000000000001 2 8 func ts s 1 3951 s simpson 0 0000000000000000001 2 16 func ts s 1 4159 s simpson 0 0000000000000000001 2 32 func ts s 1 4263 所以可以得到 S4 1 3534 S8 1 3951 S16 1 4159 S32 1 4263 3 待解方程 function y func ts x y 1 1 sin x 2 运行结果 s simpson 0 pi 2 4 func ts s 1 1107 s simpson 0 pi 2 8 func ts 11 s 1 1107 s simpson 0 pi 2 16 func ts s 1 1107 s simpson 0 pi 2 32 func ts s 1 1107 所以可以得到 S4 1 1107 S8 1 1107 S16 1 1107 S32 1 1107 五 常微分方程五 常微分方程 1 试用 Euler 法及改进 Euler 法计算初值问题 1 Euler 法 待解方程 function f descent t y f y 2 t y 运行过程 clear clc a 0 b 1 h 0 2 n b a h t 1 0 y 1 1 for i 2 n 1 t i t 1 i 1 h y i y i 1 h descent t i 1 y i 1 end for i 1 n 1 fprintf t f t i fprintf y f n y i end plot t y xlabel t ylabel y 运行结果 a 0 b 1 t 0 000000 y 1 000000 t 0 200000 y 1 200000 t 0 400000 y 1 373333 t 0 600000 y 1 531495 t 0 800000 y 1 681085 t 1 000000 y 1 826948 图形 12 2 改进 Euler 法 函数程序 function f algeb x y f y 2 x y 求解过程 clear clc a 0 b 1 h 0 2 n b a h t 1 0 y 1 1 for i 2 n 1 t i t 1 i 1 h yp i y i 1 h algeb t i 1 y i 1 y i y i 1 0 5 h algeb t i 1 y i 1 algeb t i yp i end for i 1 n 1 fprintf t f t i fprintf y f n y