MATLAB源代码.doc

埃特金插值埃特金插值 function f Atken x y x0 syms t if length x length y n length x else disp return end y1 1 n t for i 1 n 1 for j i 1 n y1 j y j t x i x j x i y i t x j x i x j end y y1 simplify y1 end if nargin 3 f subs y1 n t x0 else simplify y1 n f collect y1 n f vpa f 6 end 特特征征多多项项式式法法 function l Chapoly A 特征多项式法求矩阵特征值 已知矩阵 A 求得的矩阵特征值 l syms t N size A n N 1 1 y det A t eye n n l solve y l vpa 1 5 function f Chebyshev y k x0 syms t T 1 k 1 t T 1 1 T 2 t c 1 k 1 0 0 c 1 int subs y findsym sym y sym t T 1 sqrt 1 t 2 t 1 1 pi c 2 2 int subs y findsym sym y sym t T 2 sqrt 1 t 2 t 1 1 pi f c 1 c 2 t for i 3 k 1 T i 2 t T i 1 T i 2 c i 2 int subs y findsym sym y sym t T i sqrt 1 t 2 t 1 1 2 f f c i T i f vpa f 6 if i k 1 if nargin 3 f subs f t x0 else f vpa f 6 end end end 弦割法求解非线性方程弦割法求解非线性方程 function x k t ChordsecantToEquation f x0 x1 eps 弦割法求解非线性方程 x k t ChordsecantToEquation f x0 x1 eps x 近似解 k 迭代次数 t 运算时间 f 原函数 定义为内联函数 x0 x1 初始值 eps 误差限 应用举例 f inline x 3 4 x 2 10 x ChordsecantToEquation f 1 2 0 5e 6 x k ChordsecantToEquation f 1 2 0 5e 6 x k t ChordsecantToEquation f 1 2 0 5e 6 函数的最后一个参数也可以不写 默认情况下 eps 0 5e 6 x k t ChordsecantToEquation f 1 2 if nargin 3 eps 0 5e 6 end tic k 0 while 1 x x1 f x1 x1 x0 f x1 f x0 k k 1 if abs x x1 30 break end x0 x1 x1 x end t toc if k 30 disp 迭代次数太多 x 0 t 0 end 复合求积公式复合求积公式 function I step CombineTraprl f a b eps f 被积函数 a b 积分上下限 eps 精度 I 积分结果 step 积分的子区间数 if nargin 3 eps 1 0e 4 end n 1 h b a 2 I1 0 I2 subs sym f findsym sym f a subs sym f findsym sym f b h while abs I2 I1 eps n n 1 h b a n I1 I2 I2 0 for i 0 n 1 x a h i x1 x h I2 I2 h 2 subs sym f findsym sym f x subs sym f findsym sym f x1 end end I I2 step n function f DCS x y x0 syms t if length x length y n length x c 1 n 0 0 else disp x 和 y 的维数不相等 return end c 1 y 1 for i 1 n 1 for j i 1 n y1 j x j x i y j y i end c i 1 y1 i 1 y y1 end f c n for i 1 n 1 f c n i t x n i f f vpa f 6 if i n 1 if nargin 3 f subs f t x0 else f vpa f 6 end end end 用欧拉法求一阶常微分方程的数值解用欧拉法求一阶常微分方程的数值解 function y DEEuler f h a b y0 varvec format long N b a h y zeros N 1 1 x a h b y 1 y0 for i 2 N 1 y i y i 1 h Funval f varvec x i 1 y i 1 end format short 用隐式欧拉法求一阶常微分方程的数值解用隐式欧拉法求一阶常微分方程的数值解 function y DEimpEuler f h a b y0 varvec format long N b a h y zeros N 1 1 y 1 y0 x a h b var findsym f for i 2 N 1 fx Funval f var 1 x i gx y i 1 h fx varvec 2 y i NewtonRoot gx 10 10 eps end format short 龙格库塔法龙格库塔法 function y DELGKT2 suen f h a b y0 varvec format long N uint16 b a h y zeros N 1 1 y 1 y0 x a h b var findsym f for i 2 N 1 K1 Funval f varvec x i 1 y i 1 K2 Funval f varvec x i 1 2 h 3 y i 1 K1 2 h 3 y i y i 1 h K1 3 K2 4 end format short 用改进欧拉法求一阶常微分方程的数值解用改进欧拉法求一阶常微分方程的数值解 function y DEModifEuler f h a b y0 varvec format long N b a h y zeros N 1 1 y 1 y0 x a h b var findsym f for i 2 N 1 v1 Funval f varvec x i 1 y i 1 t y i 1 h v1 v2 Funval f varvec x i t y i y i 1 h v1 v2 2 end format short 一阶微分方程一阶微分方程 Function dy myfun 1 x y 定义输入 输出变量和函数文件名 dy zeros 2 1 明确 dy 的维数 用微分方程组时不可缺省 dy 1 y 2 dy m 表示 y 的 m 阶导数 y n 表示 y 的第 n 列 dy 2 2 y 2 2 y 1 5 exp 2 x sin x 与方程组中第二个微分方程对应 基本基本 GuassGuass 消去法求解线性方程组消去法求解线性方程组 function x EqtsBasicGuass A b 基本 Guass 消去法求解线性方程组 Ax b x EqtsBasicGuass A b x 解向量 列向量 A 线性方程组的矩阵 b 列向量 检查输入参数 if size A 1 size b 1 disp 输入参数有误 x return end A b A A b 消去过程 n size A 1 l zeros n for k 1 n 1 for i k 1 n l i k A i k A k k end for i k 1 n for j k 1 n 1 A i j A i j l i k A k j end for j 1 k A i j 0 end end end 回代过程 x zeros n 1 x n A n n 1 A n n for i n 1 1 1 y 0 for j i 1 n y y A i j x j end x i A i n 1 y A i i end return 追赶法求解三对角线性方程组追赶法求解三对角线性方程组 function x EqtsForwardAndBackward L D U b 追赶法求解三对角线性方程组 Ax b x EqtsForwardAndBackward L D U b x 三对角线性方程组的解 L 三对角矩阵的下对角线 行向量 D 三对角矩阵的对角线 行向量 U 三对角矩阵的上对角线 行向量 b 线性方程组 Ax b 中的 b 列向量 应用举例 L 1 2 3 D 2 3 4 5 U 1 2 3 b 6 1 2 1 x EqtsForwardAndBackward L D U b 检查参数的输入是否正确 n length D m length b n1 length L n2 length U if n n1 1 n n2 1 n m disp 输入参数有误 x return end 追的过程 for i 2 n L i 1 L i 1 D i 1 D i D i L i 1 U i 1 end x zeros n 1 x 1 b 1 for i 2 n x i b i L i 1 x i 1 end 赶的过程 x n x n D n for i n 1 1 1 x i x i U i x i 1 D i end return Guass Seidel 迭代法迭代法 function x k EqtsGS A b x0 eps Guass Seidel 迭代法求解线性方程组 Ax b x k EqtsGS A b x0 eps x 解向量 列向量 k 迭代次数 A 系数矩阵 b 列向量 x0 迭代初始值 列向量 eps 误差限 可缺省 缺省值为 0 5e 6 应用举例 A 10 3 1 2 10 3 1 3 10 b 14 5 14 x0 0 0 0 x k EqtsGS A b x0 0 5e 6 x EqtsGS A b x0 if nargin 3 eps 0 5e 6 end 检查输入参数 n length b if size A 1 n n length x0 disp 输入参数有误 x k return end 迭代求解 k 0 x zeros n 1 while 1 k k 1 for i 1 n z 0 for j 1 i 1 z z A i j x j end for j i 1 n z z A i j x0 j end x i b i z A i i end if norm x x0 30 break end x0 x end if k 30 disp 迭代次数太多 x return end return JacobiJacobi 迭代法求解线性方程组迭代法求解线性方程组 function x k EqtsJacobi A b x0 eps Jacobi 迭代法求解线性方程组 Ax b x k EqtsJacobi A b x0 eps x 解向量 列向量 k 迭代次数 A 系数矩阵 b 列向量 x0 迭代初始值 列向量 eps 误差限 可缺省 缺省值为 0 5e 6 应用举例 A 10 3 1 2 10 3 1 3 10 b 14 5 14 x0 0 0 0 x k EqtsJacobi A b x0 0 5e 6 x EqtsJacobi A b x0 if nargin 3 eps 0 5e 6 end 检查输入参数 n length b if size A 1 n n length x0 disp 输入参数有误 x k return end 迭代求解 k 0 x zeros n 1 while 1 k k 1 for i 1 n z 0 for j 1 i 1 z z A i j x0 j end for j i 1 n z z A i j x0 j end x i b i z A i i end if norm x x0 30 break end x0 x end if k 30 disp 迭代次数太多 x return end return 逐次超松驰迭代法逐次超松驰迭代法 function x k EqtsSOR A b x0 omiga eps 超松弛 SOR Successive Over Relaxation 迭代法求解线性方程组 Ax b x k EqtsSOR A b x0 eps x 解向量 列向量 k 迭代次数 A 系数矩阵 b 列向量 x0 迭代初始值 列向量 omiga 松弛因子 可缺省 缺省值为 1 即为 GS 迭代法 eps 误差限 可缺省 缺省值为 0 5e 6 应用举例 A 4 3 0 3 4 1 0 1 4 b 24 30 24 x0 1 1 1 omiga 1 25 x k EqtsSOR A b x0 omiga 0 5e 6 x EqtsSOR A b x0 if nargin 4 eps 0 5e 6 end if nargin 3 omiga 1 eps 0 5e 6 end 检查输入参数 n length b if size A 1 n n length x0 disp 输入参数有误 x k return end 迭代求解 k 0 x zeros n 1 while 1 k k 1 for i 1 n z 0 for j 1 i 1 z z A i j x j end for j i 1 n z z A i j x0 j end x i 1 omiga x0 i omiga b i z A i i end if norm x x0 30 break end x0 x end if k 30 disp 迭代次数太多 x return end return 连分式插值连分式插值 function f DCS x y x0 syms t if length x length y n length x c 1 n 0 0 else disp return end c 1 y 1 for i 1 n 1 for j 1 1 n y1 j x j x i y j y i end c i 1 y1 i 1 y y1 end f c n for i 1 n 1 f c n i t x n i f f vpa f 6 if i n 1 if nargin 3 f subs f t x0 else f vpa f 6 end end end 用二分法求方程的一个根用二分法求方程的一个根 function root HalfInterval f a b eps if nargin 3 eps 1 0e 4 end f1 subs sym f findsym sym f a f2 subs sym f findsym sym f b if f1 0 root a end if f2 0 root b end if f1 f2 0 disp 两端点函数值乘积大于 0 return else root FindRoots f a b eps end function r FindRoots f a b eps f 1 subs sym f findsym sym f a f 2 subs sym f findsym sym f b mf subs sym f findsym sym f a b 2 if f 1 mf 0 t a b 2 r FindRoots f t b eps else if f 1 mf 0 r a b 2 else if abs b a eps r b 3 a 4 else s a b 2 r FindRoots f a s eps end end end HermiteHermite 插值的插值的 matlabmatlab 代码代码 function f Hermite x y y 1 x0 syms t f 0 0 if length x length y if length y length y 1 n length x else disp y 和 y 的导数的维数不相等 return end else disp x 和 y 的维数不相等 return end for i 1 n h 1 0 a 0 0 for j 1 n if j i h h t x j 2 x i x j 2 a a 1 x i x j end end f f h x i t 2 a y i y 1 i y i if i n if nargin 4 f subs f t x0 else f vpa f 6 end end end 辛普生公式的求积辛普生公式的求积 function S n IntCompSimpson f D4f a b eps if nargin 4 eps 0 5e 7 默认精度 end 求被积函数的四次导数在区间 a b 上的最大值 x fval fminbnd D4f a b optimset TolX eps 10 fmax fval 计算等分区间数 n ceil sqrt sqrt abs b a 5 fmax 16 180 eps h b a n 步长 S f a f b 4 f a h 2 for k 1 n 1 x1 a k h x2 x1 h 2 S S 4 f x2 2 f x1 end S S h 6 return 复化梯形公式求积分复化梯形公式求积分 function T n IntCompTrape f D2f a b eps 复化梯形公式求积分 T n IntCompTrape f D2f a b eps T 求积结果 n 区间等分数 f 被积函数 可利用函数脚本文件事先定义 也可以利用内联函数 f 被积函数的二次导数的相反值 可利用函数脚本文件事先定义 也可以利用内联函数 取相反值是为了便于计算被积函数的二次导数在区间 a b 上的最大值 a 积分下限 b 积分上限 eps 误差限 应用举例 事先定义 T n IntCompTrape IntF IntD2F 0 1 T n IntCompTrape IntF IntD2F 0 1 0 5e 7 利用内联函数 F inline 3 x D2F inline log 3 2 3 x T n IntCompTrape F D2F 0 1 if nargin 4 eps 0 5e 7 默认精度 end 求被积函数的二次导数在区间 a b 上的最大值 x fval fminbnd D2f a b optimset TolX eps 10 fmax fval 计算等分区间数 n ceil sqrt abs fmax b a 3 12 eps h b a n 步长 T f a f b for k 1 n 1 x1 a k h T T 2 f x1 end T h T 2 return 拉格朗日插值拉格朗日插值 1 function f InterpLagrange2 x y x0 构造 Lagrange 插值多项式 此函数中借助符号运算来求 Lagrange 基函数 运算速度较慢 不推荐此种方法 f InterpLagrange2 x y x0 f 插值多项式或者是插值多项式在 x0 处的值 x 节点 y 函数值 x0 某一测试点 调用格式 f InterpLagrange2 x y 返回插值多项式 f InterpLagrange2 x y x0 返回插值多项式在点 x0 处的值 举例 x 0 32 0 34 0 36 y 0 314567 0 333487 0 352274 x0 0 33 f InterpLagrange2 x y f InterpLagrange2 x y x0 if length x length y n length x else disp 节点个数和函数值个数不同 f return end syms t f 0 for i 1 n l y i for j 1 n if j i continue end l l t x j x i x j 借助符号运算 计算 Lagrange 基函数 end f f l simplify f 化简多项式 if i n if nargin 3 f subs f t x0 计算插值多项式 f 在点 x0 处的值 else f collect f 计算插值多项式 展开并合并同类项 f vpa f 6 设置多项式系数的有效数字 end end end 拉格朗日插值拉格朗日插值 2 function f Language x y x0 syms t if length x length y n length x else disp return end f 0 0 for i 1 n l y i for j 1 i 1 l l t x j x i x j end for j i 1 n l l t x j x i x j end f f 1 simplify f if i n if nargin 3 f subs f t x0 else f collect f f vpa f 6 end end end 勒让德多项式勒让德多项式 function f Legendre y k x0 syms t P 1 k 1 t P 1 1 P 2 t c 1 k 1 0 0 c 1 int subs y findsym sym y sym t P 1 t 1 1 2 c 2 int subs y findsym sym y sym t P 2 t 1 1 2 f c 1 c 2 t for i 3 k 1 P i 2 i 3 P i 1 t i 2 P i 2 i 1 c i int subs y findsym sym y t P i t 1 1 2 f f c i P i if i k 1 if nargin 3 f subs f t x0 else f vpa f 6 end end end 高斯高斯 勒让德求积公式勒让德求积公式 function I IntGauss f a b n AK XK if n 5 XK 0 else XK1 b a 2 XK a b 2 I b a 2 sum AK subs sym f findsym f XK1 end ta b a 2 tb a b 2 switch n case 0 I 2 ta subs sym f findsym sym f tb case 1 I ta subs sym f findsym sym f ta 0 5773503 tb subs sym f findsym sym f ta 0 5773503 tb case 2 I ta 0 55555556 subs sym f findsym sym f ta 0 7745967 tb 0 55555556 subs sym f findsym sym f ta 0 7745967 tb 0 88888889 subs sym f findsym sym f tb case 3 I ta 0 3478548 subs sym f findsym sym f ta 0 8611363 tb 0 3478548 subs sym f findsym sym f ta 0 8611363 tb 0 6521452 subs sym f findsym sym f ta 0 3398810 tb 0 6521452 subs sym f findsym sym f ta 0 3398810 tb case 4 I ta 0 2369269 subs sym f findsym sym f ta 0 9061793 tb 0 2369269 subs sym f findsym sym f ta 0 9061793 tb 0 4786287 subs sym f findsym sym f ta 0 5384693 tb 0 4786287 subs sym f findsym sym f ta 0 5384693 tb 0 5688889 subs sym f findsym sym f tb End LU 分解 分解 1 function L U flag LU Decom A n m size A 要求所分解的矩阵是方阵 否则停止计算 并输出出错信息 if n m error The rows and columns of matrix A must be equal return end L eye n U zeros n flag OK for k 1 n for j k n z 0 for q 1 k 1 z z L k q U q j end U k j A k j z end if abs U k k dfb 初始值取两端点导数较大者 root a fa dfa else root b fb dfb end while tol eps r1 root fx subs sym f findsym sym f r1 dfx subs sym fun findsym sym fun r1 求该点的导数值 root r1 fx dfx 迭代的核心公式 tol abs root r1 end 用抛物线法求方程的一个根用抛物线法求方程的一个根 function root Parabola f a b x eps if nargin 4 eps 1 0e 4 end f1 subs sym f findsym sym f a f2 subs sym f findsym sym f b if f1 0 root a end if f2 0 root b end if f1 f2 0 disp 两端点函数值乘积大于 0 return else tol 1 fa subs sym f findsym sym f a fb subs sym f findsym sym f b fx subs sym f findsym sym f x d1 fb fa b a d2 fx fb x b d3 d2 d1 x a B d2 d3 x b root x 2 fx B sign B sqrt B 2 4 fx d3 t zeros 3 t 1 a t 2 b t 3 x while tol eps t 1 t 2 t 2 t 3 t 3 root f1 subs sym f findsym sym f t 1 f2 subs sym f findsym sym f t 2 f3 subs sym f findsym sym f t 3 d1 f2 f1 t 2 t 1 d2 f3 f2 t 3 t 2 d3 d2 d1 t 3 t 1 B d2 d3 t 3 t 2 root t 3 2 f3 B sign B sqrt B 2 4 f3 d3 tol abs root t 3 end end 幂法求矩阵的主特征值及主特征向量幂法求矩阵的主特征值及主特征向量 function l v s pmethod A x0 eps if nargin 2 eps 1 0e 6 end v x0 v 为主特征向量 M 5000 迭代步数限制 m 0 l 0 for k 1 M y A v m max y m 为按模最大的分量 v y m if abs m l eps l m 到所需精度 退出 l 为主特征值 s k s 为迭代步数 return else if k M disp 迭代步数太多 收敛速度太慢 l m s M else l m end end end ipmethod 收缩法求矩阵全部特征值收缩法求矩阵全部特征值 function l v ipmethod A x0 eps if nargin 2 eps 1 0e 6 end v x0 M 5000 迭代步数限制 m 0 l 0 for k 1 M y A v m max y m 为按模最大的分量 v y m if abs m l eps l 1 m 到所需精度 退出 l 为按模的最小特征值 return else if k M disp 迭代步数太多 收敛速度太慢 l 1 m else l m end end end QR 基本算法求矩阵全部特征值基本算法求矩阵全部特征值 function l qrtz A M for i 1 M M 为迭代步数 q r qr A A r q l diag A end 龙格贝格求积公式龙格贝格求积公式 function T romberg a b detal a 是积分下限 b 是积分上限 detal 是要 求误差 n 1 n 的初值为 1 while 1 T zeros n 1 n 1 构造 n 1 n 1 的矩阵 T T 1 1 b a feval f a feval f b 2 f 为被积函数 for k 1 n h b a 2 k m 2 k 1 s 0 for i 1 m s s feval f a h 2 i 1 end T k 1 1 T k 1 2 h s 算出矩阵 T 的第一列 end for p 1 n for j p n T j 1 p 1 4 p T j 1 p T j p 4 p