MATLAB上机作业(数值分析方法(上理))

解得 x=0.8776 牛顿法求解：
建立第二个函数文件 function f2=f2(x) f2=15*x^2+12*x;
计算代码
x(1)=-1.7;jingdu=0.002; x(2)=x(1)-(f(x(1))/f2(x(1)));n=2;
while abs(x(n)-x(n-1))>jingdu x(n+1)=x(n)-(f(x(n))/f2(x(n))); n=n+1; end x
自定义计算文件
X = [0 pi/6 pi/4 pi/3 pi/2]; Y = [0 0.5 0.7071 0.8660 1]; x = linspace(0,pi,50); M = 1; [y, R] = lagrange(X, Y, x, M); y1 = sin(x); %画图 errorbar(x,y,R,'.k') hold on plot(X, Y, 'or', x, y, '*r', x, y1, '-b'); axis([-0.5,3.5,-1,1.2]);%把坐标边界加大一点，不然坐标点遮 住坐标轴了
end
A(i)=x^(i-1);
% 计算 G for j=1:n+1 for i=1:n+1 for k=1:m+1 B(k)=W(k)*subs(A(i),x,k)*subs(A(j),x,k); end G(i,j)=sum(B); end end G % 计算 d for j=1:n+1 for k=1:m+1 B(k)=W(k)*F(k)*subs(A(j),x,k); end d(j)=sum(B); end d=d' % 求出拟合曲线的系数 G=G^-1; C=G*d; S=0; for i=1:n+1 S=S+C(i)*x^(i-1); end % 画图 scatter(X,F,'*b'); hold on; ezplot(S,[0,6]); xlabel('x') ylabel('y') title('拟合曲线') end
end C = A(n, n); q1 = abs(q1*(z-X(n))); for k = (n-1):-1:1 C = conv(C, poly(X(k))); d = length(C); C(d) = C(d) + A(k,k);%在最后一维，也就是常数项加 上新的差商 end y(t) = polyval(C,z);%向量 y 是向量 x 处的插值，误差限 R，n 次牛顿插值多项式 L 及其系数向量 C end L = poly2sym(C); R(t) = M * q1 / c1;
上述代码插值部分和数值积分均参考自 CSDN 的几位知名博主， 大家可以自行查找。本文档旨在分享我的数值分析上机作业，不 做商用，故不与人产生任商业用途，转载请勿商用。想要系统 学习 matlab 软件，推荐飞天小苗苗的博客以及优酷她的教学视 频。
《数值计算方法》 上机报告
班级：机械*班 学号：171441**** 姓名：范涛
目 录 1. 非线性方程的数值解法 1.1 问题描述 1.2 基本原理和步骤 1.3 程序源代码 1.4 计算结果及分析 1.5 小结 2. 函数插值与曲线拟合 2.1 问题描述 2.2 基本原理和步骤 2.3 程序源代码 2.4 计算结果及分析 2.5 小结 3. 数值积分 3.1 问题描述 3.2 基本原理和步骤 3.3 程序源代码 3.4 计算结果及分析 3.5 小结
解得 x=0.8776 割线法求解： 计算代码
x(1)=-1.7;jingdu=0.002; x(2)=x(1)-(f(x(1))/f2(x(1)));n=2; while abs(x(n)-x(n-1))>jingdu x(n+1)=x(n)-(f(x(n))/(f(x(n))-f(x(n-1))))*(x(n)-x( n-1)); n=n+1; end x
解得 x=0.8776 1.4 计算结果分析 各种方法只要能够求解的情况下，多次求解就可以无限逼近所求解， 但是迭代法必须要求在一定范围内有极限， 而割线法要求一定范围内 必须为可切曲线，故有局限性 2.函数插值与曲线拟合 2.1 问题描述
用 Lagrange 插值来求 sinx 在某点的值，并估计其误差，已知 sin0° = 0, sin30° = 0.5, sin45° = 0.7071, sin60° = 0.8660, sin90° = 1.
解得 ans =18.849555921538759 复化三点式（辛普森公式）
function Sn = Sn(a,b,n) format long h = (b-a)/n; sum1 = 0; sum2 = 0; for i = 0:n-1 sum1 = sum1 + f(a+(i+1/2).*h); end for j = 1:n-1 sum2 = sum2 + f(a+j.*h); end Sn = h/6*(f(a)+4*sum1+2*sum2+f(b));
割线法：是利用牛顿迭代法的思想,在根的某个领域内，函数有直至 二阶的连续导数，并且不等于 0，则在领域内选取初值 x0,x1,迭代 均收敛。 1.3 程序源代码 建立被求解函数文件 function f=f(x) f=5*x^3+6*x^2-8; 二分法计算代码： a=-2;b=1;jingdu=0.002 while abs(a-b)>jingdu x=(a+b)/2;
c1 = c1 * j;
自定义计算文件
X = [0 pi/6 pi/4 pi/3 pi/2];%X 是 n+1 个节点(x_i,y_i)(i = 1,2, ... , n+1)横坐标向量，Y 是纵坐标向量 Y = [0 0.5 0.7071 0.8660 1]; x = linspace(0,pi,50); %x 是以向量形式输入的 m 个插值点，M
1.非线性方程的数值解法 1.1 令 f(x)=5*x^3+6*x^2-8=0，求根,x∈[-2,1],ε=0.002 1.2 基本原理和步骤 二分法：将[a,b]分为相等的 2 个小区间，计算小区间端点函数值， 找出新的有根区间。重复上述过程，直到找到满足精度要求的根。 迭代法：是取 x0 之后，在这个基础上，找到比 x0 更接近的方程的 跟，一步一步迭代，从而找到更接近方程根的近似跟。 牛顿法： 把非线性函数 f(x)=0 在 x0 处展开成泰勒级数取其线性 部分，作为非线性方程的近似方程， 则有 设 ，则其解为
legend('误差','样本点','拉格朗日插值估算','sin(x)'); title('拉格朗日插值法拟合函数');
牛顿插值法： 定义牛顿差值多项式函数
function[y,R,A,C,L] = newton(X,Y,x,M) n = length(X); m = length(x); for t = 1 : m z = x(t); A = zeros(n,n);%差商的矩阵 A A(:,1) = Y'; s = 0.0; p = 1.0; q1 = 1.0; c1 = 1.0; for j = 2 : n for i = j : n A(i,j) = (A(i,j-1) A(i-1,j-1))/(X(i)-X(i-j+1)); end q1 = abs(q1*(z-X(j-1)));
解得 ans =18.849555921538759 3.4 计算结果及分析 原函数需要多次积分能得到不定积分，应该是我的五次方对于我 的区间选择的影响较大，使得五次方之后最终结果前面的有效数 字接近，尝试把要求的函数降次积分之后返现结果误差较大，是 题目的问题，程序正确。 复化 Cotes 公式的代数精度最高， 复化 Simpson 其次， 复化梯形 公式的代数精度最低
2.2 基本原理和步骤 2.3 程序源代码 拉格朗日插值法： 拉格朗日插值函数文件
function [y, R] = lagrange(X, Y, x, M) n = length(X);m = length(x); for i = 1:m z = x(i); s = 0.0; for k = 1:n p = 1.0; q1 = 1.0; c1 = 1.0; for j = 1:n if j~=k p = p * (z - X(j)) / (X(k) - X(j)); end q1 = abs(q1 * (z - X(j))); c1 = c1 * j; end s = p * Y(k) + s; end y(i) = s; R(i) = M * q1 / c1; end
解得 ans =18.849555921538766 复化五点式（柯特斯公式）
function Cn = Cn(a,b,n) format long h = (b-a)/n;
sum1 = 0; sum2 = 0; for i = 0:n-1 sum1 = sum1 + 32*f(a+(i+1/4).*h)+12*f(a+(i+1/2).*h)+32*f(a+(i+3/4). *h); end for j = 1:n-1 sum2 = sum2 + 14*f(a+j.*h); end Cn = h/90*(7*f(a)+sum1+sum2+7*f(b));
命令行输入： X = [0 pi/6 pi/4 pi/3 pi/2];
F = [0 0.5 0.7071 0.8660 1]; W=[1 1 1 1 1]; S=mypolyfit(X,F,W,4,3)