牛顿插值法matlab代码

合集下载

数值计算实验牛顿前插和后插插值

数值计算实验牛顿前插和后插插值

数值计算实验二姓名:方小开学号:20060810202 班级:计科0602一. 实验目的:1、差分的matlab实现;2、Newton插值的matlab实现;二. 实验原理:MATLAB在线性代数,矩阵分析,数值及优化,数理统计和随机信号分析,电路系统,系统动力学,信号与图像处理,控制理论分析和系统设计,过程控制,建模和仿真,通信系统,等有广泛的应用。

它具有功能强大,界面友好,语言自然即开放性等特点。

三.试验环境MATLAB7.0四. 试验过程及现象:1、牛顿插值公式:把下面的matlab程序在matlab中建立M-file文件并保存;function [d,v1]=newtonjz(x,y,v) %d 插商表 v是要插入x v1是插入的y值n=length(x);d=zeros(n,n);d(:,1)=y';for j=2:nfor i=j:nd(i,j)=(d(i,j-1)-d(i-1,j-1))/(x(i)-x(i-j+1));endendw=1;v1=d(1,1);for i=2:nw=w*(v-x(i-1));v1=v1+d(i,i)*w;end分别给x,y赋初值,并调用Newton插值函数得到结果如下:x=[0.40,0.55,0.65,0.80,0.90,1.05];y=[0.41075,0.57815,0.69675,0.88811,1.02652,1.25382];[z,xy]=newtonjz(x,y,0.596);z =0.4108 0 0 0 0 00.5782 1.1160 0 0 0 00.6967 1.1860 0.2800 0 0 00.8881 1.2757 0.3589 0.1973 0 01.0265 1.3841 0.4335 0.2130 0.0312 01.2538 1.5153 0.5249 0.2287 0.0314 0.0003xy =0.63192Newton前插公式:把Newton前插公式的matlab程序写在matlab中建立M-file文件并保存;function [d,v1]=newtonBefore(x,y,t)n=length(x);d=zeros(n,n);d(:,1)=y';for j=2:nfor i=1:n-j+1d(i,j)=(d(i+1,j-1)-d(i,j-1));endendw=1;m=1;v1=d(1,1);for i=2:nw=w*(t-i+2);m=m*(i-1);v1=v1+d(1,i)*(w/m);end分别给x,y赋初值,并调用Newton前插函数得到结果如下;x=[1 1.05 1.10 1.15 1.20 1.25 1.30];y=[1 1.0247 1.04881 1.07238 1.09544 1.11803 1.14017];>> [z,qc]=newtonBefore(x,y,0.2);qc =1.004992263808003、Newton后插公式:把Newton后插公式的matlab程序写在matlab中建立M-file文件并保存;function [d,v1]=newtonAfter(x,y,t)n=length(x);d=zeros(n,n);d(:,1)=y';for j=2:nfor i=j:nd(i,j)=(d(i,j-1)-d(i-1,j-1));endendw=1;m=1;v1=d(1,1);for i=2:nw=w*(t+i-2);m=m*(i-1);v1=v1+d(n,i)*(w/m);end分别给x,y赋初值,并调用Newton后插函数得到结果如下;x=[1 1.05 1.10 1.15 1.20 1.25 1.30];y=[1 1.0247 1.04881 1.07238 1.09544 1.11803 1.14017];>> [z,hc]=newtonAfter(x,y,-0.4);hc=1.13136982835200五.遇到的问题在调试的过程中也遇到了一些小小的问题,如输出的结果只显示4位有效数字,结果的精度太低了,不能满足要求,因此在matlab中把数据的格式从short型改成了long型,这样就大大的提高了结果的精确度,减少了误差。

牛顿-柯特斯公式matlab

牛顿-柯特斯公式matlab

牛顿-柯特斯公式matlab首先,让我们来了解一下数值积分的基本概念。

数值积分是通过求取一个函数在给定区间上的近似面积来计算函数的定积分。

一种常见的数值积分方法是使用插值多项式来近似函数,并在给定区间上对该多项式进行积分。

牛顿插值多项式是由一组不同的x值和对应的函数值构成的。

该多项式通过这些点来逼近函数,并可以用于在任意点上计算函数的近似值。

牛顿插值多项式的形式如下:P(x)=f[x₀]+f[x₀,x₁](x-x₀)+f[x₀,x₁,x₂](x-x₀)(x-x₁)+...其中,f[x₀]表示函数在x₀上的值,f[x₀,x₁]表示函数在x₀和x₁上的差商。

柯特斯系数用于计算牛顿插值多项式在给定区间上的积分。

公式如下:C₀=1C₁=h/2C₂=h²/6C₃=h³/12C₄=h⁴/20其中,h表示区间的宽度。

在MATLAB中,可以使用以下代码来实现牛顿-柯特斯公式:function result = newton_cotes(f, a, b, n)h=(b-a)/n;x=a:h:b;fx = f(x);coefficient = zeros(n+1, 1);coefficient(1) = 1;for i = 2:n+1coefficient(i) = coefficient(i-1) * (h^(i-1)) / factorial(i-1);endresult = sum(fx .* coefficient);end```在上面的代码中,`f`表示要积分的函数,`a`和`b`表示积分区间的起始点和结束点,`n`表示节点的数量。

首先,我们计算出节点的横坐标和对应的函数值。

然后,根据柯特斯系数的公式计算系数。

最后,将函数值与系数相乘,并求和,从而得到近似的积分值。

例如,我们要计算函数f(x) = sin(x)在区间[0, π/2]上的积分值,可以使用以下代码:a=0;b = pi/2;n=4;result = newton_cotes(f, a, b, n);disp(result);```运行该代码,将输出函数f(x)在区间[0,π/2]上的近似积分值。

拉格朗日插值和牛顿插值matlab

拉格朗日插值和牛顿插值matlab

拉格朗⽇插值和⽜顿插值matlab1. 已知函数在下列各点的值为0.20.40.60.8 1.00.980.920.810.640.38⽤插值法对数据进⾏拟合,要求给出Lagrange插值多项式和Newton插值多项式的表达式,并计算插值多项式在点的值。

程序:x=[0.2 0.4 0.6 0.8 1.0];y=[0.98 0.92 0.81 0.64 0.38];x0=[0.2 0.28 0.44 0.76 1 1.08];[f,f0]=Lagrange(x,y,x0)function [f,f0] = Lagrange(x,y,x0)%求已知数据点的Lagrange插值多项式f,并计算插值多项式f在数据点x0的函数值f0syms t;n = length(x);f = 0.0;for i = 1:nl = y(i);for j = 1:i-1l = l*(t-x(j))/(x(i)-x(j));end;for j = i+1:nl = l*(t-x(j))/(x(i)-x(j));end;f = f + l;simplify(f);if(i==n)f0 = subs(f,'t',x0);f = collect(f);f = vpa(f,6);endend结果:>> Untitled3f =- 0.520833*t^4 + 0.833333*t^3 - 1.10417*t^2 + 0.191667*t + 0.98f0 =[ 49/50, 60137/62500, 56377/62500, 42497/62500, 19/50, 15017/62500]⽜顿:%y为对应x的值,A为差商表,C为多项式系数,L为多项式%X为给定节点,Y为节点值,x为待求节点function[y,A,C,L] = newton(X,Y,x,M)n = length(X);m = length(x);for t = 1 : mz = x(t);A = zeros(n,n);A(:,1) = Y';s = 0.0; p = 1.0; q1 = 1.0; c1 = 1.0;for j = 2 : nfor i = j : nA(i,j) = (A(i,j-1) - A(i-1,j-1))/(X(i)-X(i-j+1));endq1 = abs(q1*(z-X(j-1)));c1 = c1 * j;endC = A(n, n); q1 = abs(q1*(z-X(n)));for k = (n-1):-1:1C = conv(C, poly(X(k)));d = length(C);C(d) = C(d) + A(k,k);endy(t) = polyval(C,z);endL = poly2sym(C);x=[0.2 0.4 0.6 0.8 1.0];y=[0.98 0.92 0.81 0.64 0.38];x0=[0.2 0.28 0.44 0.76 1 1.08];m=1;[y,A,C,L]=newton(x,y,x0,m)结果:y =0.9800 0.9622 0.9020 0.6800 0.3800 0.2403A =0.9800 0 0 0 00.9200 -0.3000 0 0 00.8100 -0.5500 -0.6250 0 00.6400 -0.8500 -0.7500 -0.2083 00.3800 -1.3000 -1.1250 -0.6250 -0.5208C =-0.5208 0.8333 -1.1042 0.1917 0.9800L =- (25*x^4)/48 + (5*x^3)/6 - (53*x^2)/48 + (23*x)/120 + 49/502. 在区间上分别取,⽤两组等距节点对Runge函数作多项式插值(Lagrange插值和Newton插值均可),要求对每个值,分别画出插值多项式和函数的曲线。

matlab编程实现二分法牛顿法黄金分割法最速下降matlab程序代码

matlab编程实现二分法牛顿法黄金分割法最速下降matlab程序代码

matlab编程实现二分法牛顿法黄金分割法最速下降matlab程序代码二分法(Bisection Method)是一种寻找函数零点的数值计算方法。

该方法的基本思想是:首先确定一个区间[a, b],使得函数在这个区间的两个端点处的函数值异号,然后将区间逐步缩小,直到找到一个区间[a', b'],使得函数在这个区间的中点处的函数值接近于零。

以下是使用MATLAB实现二分法的示例代码:```matlabfunction [x, iter] = bisection(f, a, b, tol)fa = f(a);fb = f(b);if sign(fa) == sign(fb)error('The function has the same sign at the endpoints of the interval');enditer = 0;while (b - a) / 2 > tolc=(a+b)/2;fc = f(c);if fc == 0break;endif sign(fc) == sign(fa)a=c;fa = fc;elseb=c;fb = fc;enditer = iter + 1;endx=(a+b)/2;end```牛顿法(Newton's Method)是一种用于寻找函数零点的数值计算方法。

该方法的基本思想是:通过迭代来逼近函数的零点,每次迭代通过函数的切线来确定下一个近似值,直到满足收敛条件。

以下是使用MATLAB实现牛顿法的示例代码:```matlabfunction [x, iter] = newton(f, df, x0, tol)iter = 0;while abs(f(x0)) > tolx0 = x0 - f(x0) / df(x0);iter = iter + 1;endx=x0;end```黄金分割法(Golden Section Method)是一种用于寻找函数极值点的数值计算方法。

牛顿法matlab程序及例题

牛顿法matlab程序及例题

牛顿法matlab程序及例题牛顿法是一种求解非线性方程组的常用方法,它的基本思想是通过迭代逐步逼近方程组的根。

在matlab中,可以通过编写相应的程序来实现牛顿法,并且可以通过一些例题来深入理解其应用。

下面是一份牛顿法的matlab程序:function [x, fval, exitflag, output] = mynewton(fun, x0, tol, maxiter)% fun:非线性方程组的函数句柄% x0:初始点% tol:允许误差% maxiter:最大迭代次数x = x0;fval = feval(fun, x);iter = 0;output = [];while norm(fval) > tol && iter < maxiteriter = iter + 1;J = myjacobian(fun, x);dx = - J fval;x = x + dx;fval = feval(fun, x);output = [output; [x', norm(fval)]];endif norm(fval) <= tolexitflag = 0; % 成功求解elseexitflag = 1; % 未能求解end% 计算雅可比矩阵function J = myjacobian(fun, x)n = length(x);fval = feval(fun, x);J = zeros(n);h = sqrt(eps); % 微小的增量for j = 1:nxj = x(j);x(j) = xj + h;fval1 = feval(fun, x);x(j) = xj - h;fval2 = feval(fun, x);x(j) = xj;J(:, j) = (fval1 - fval2) / (2 * h);end接下来,我们可以通过一个例题来演示牛顿法的应用。

牛顿插值法matlab程序解析

牛顿插值法matlab程序解析

牛顿插值法在MATLAB 中的实现经过n+1个不同的插值点12n+1,,x x x …,,构造牛顿插值公式1211231212n+112n =[,]()[,,]()()[,,]()()()N f x x x x f x x x x x x x f x x x x x x x x x -+--++---(x )……… 注:牛顿插值法中,用到了插值公式%我们以二次牛顿插值公式为例解析牛顿插值法的matlab 程序function[c,d ]=newpoly (x ,y )%这里x 为3个节点的横坐标组成的向量,即()123,,x x x x =,y 为纵坐标的组成向量,即()()()()123,,y f x f x f x =%c 为所得的牛顿插值多项式的系数组成的向量n=length(x);%测量向量x 的长度,即向量x 中元素i x 的个数,赋值给n,所以n=3,注:这里的“n ”仅为变量,和公式中的次数n 不一样d=zeros (n ,n ); d=zeros(3,3)%把变量d 定义为一个n 行,n 列的零矩阵,此矩阵用来储存各阶差商,格式完全等同于书中21页的表2。

1 d (:,1)=y';%此句是把向量y 的转置,即123()()()f x y f x f x ⎛⎫ ⎪= ⎪ ⎪⎝⎭,赋值给零矩阵d 的第一列%下面运用两个for 循环来构造书中21页的差商表2。

1%第一个循环(父循环),循环变量为kfor k=2:n%用来表示零矩阵d 中的第几行%第二个循环(父循环),循环变量为kfor j=k :n%用来表示零矩阵d 中的第几列d (k ,j)=(d (k ,j-1)—d(k-1,j —1))/(x (k)-x (k —j+1));%差商公式,其中d (k,j)表示零矩阵d 中的第k 行,第j 列的元素,d (k,j —1),d (k-1,j —1)等也类似,它们代表的元素随着双循环而变化,x(k —1)表示1k x -,这种计算差商的方法是根据差商表的排列位置而得来,具体解释见下面。

Langrage和Newton插值法的matlab实现

Langrage和Newton插值法的matlab实现

仅供参考1.已知数据如下:(1)用MATLAB语言编写按Langrage插值法和Newton插值法计算插值的程序,对以上数据进行插值;(2)利用MATLAB在第一个图中画出离散数据及插值函数曲线。

(1.1)langrage插值法编程实现syms xx0=[0.2,0.4,0.6,0.8,1.0];y0=[0.98,0.92,0.81,0.64,0.38];for i=1:5a=1;for j=1:5if j~=ia=expand(a*(x-x0(j)));endendb=1;for k=1:5if k~=ib=b*(x0(i)-x0(k));endendA(i)=expand(a/b);endL=0;for p=1:5L=L+y0(p)*A(p);endLL =-25/48*x^4+5/6*x^3-53/48*x^2+23/120*x+49/50(1.2)Newton插值程序实现clear allclcsyms xx0=[0.2,0.4,0.6,0.8,1.0];y0=[0.98,0.92,0.81,0.64,0.38];for k=1:5for i=1:ka=1;b=0;for j=1:kif j~=ia=a*(x0(i)-x0(j));endendb=b+y0(i)/a;endA(k)=b;endB=[1,(x-x0(1)),(x-x0(1))*(x-x0(2)),(x-x0(1))*(x-x0(2))*(x-x0(3)),(x-x 0(1))*(x-x0(2))*(x-x0(3))*(x-x0(4))];L1=A.*B;l=0;for m=1:5l=l+L1(m);endL=expand(l)L =61/100+13/30*x+383/48*x^2-155/24*x^3+475/48*x^4(2)画图x0=[0.2,0.4,0.6,0.8,1.0];y0=[0.98,0.92,0.81,0.64,0.38];subplot(1,2,1);plot(x0(1),y0(1),'+r',x0(2),y0(2),'+r',x0(3),y0(3),'+r',x0(4),y0(4),' +r',x0(5),y0(5),'+r')x=0:0.05:1;y=-25/48.*x.^4+5/6.*x.^3-53/48.*x.^2+23/120.*x+49/50;subplot(1,2,2);plot(x,y)2.给定函数21(),[1,1]125f x x x ,利用上题编好的Langrage 插值程序(或Newton 插值程序),分别取3个,5个、9个、11个等距节点作多项式插值,分别画出插值函数及原函数()f x 的图形,以验证Runge 现象、分析插值多项式的收敛性。

牛顿法matlab程序及例题

牛顿法matlab程序及例题

牛顿法matlab程序及例题牛顿法是一种求解非线性方程和优化问题的常用方法。

它利用函数的一阶和二阶导数信息来不断逼近函数的零点或极值点。

在MATLAB 中,可以用fzero函数实现非线性方程的求解,用fminunc函数实现优化问题的求解。

以下是一个简单的牛顿法的MATLAB程序示例:function [x, fx, n] = newton(f, df, x0, tol, max_iter) % f: 目标函数% df: 目标函数的一阶导数% x0: 初值% tol: 精度要求% max_iter: 最大迭代次数n = 0;while n < max_iterfx = f(x0);dfx = df(x0);if abs(dfx) < 1e-9error('牛顿法失败:一阶导数过小');endx = x0 - fx / dfx;if abs(x - x0) < tolreturn;endx0 = x;n = n + 1;enderror('牛顿法失败:达到最大迭代次数');下面是一个例题,通过牛顿法求解方程sin(x) = x / 2:f = @(x) sin(x) - x / 2;df = @(x) cos(x) - 1 / 2;[x, fx, n] = newton(f, df, 1, 1e-9, 100);fprintf('解:%.16f,函数值:%.16f,迭代次数:%d', x, fx, n);运行结果为:解:0.0000000000000000,函数值:0.0000000000000000,迭代次数:4可以看到,牛顿法很快就找到了方程的一个根。

需要注意的是,牛顿法可能会失败,特别是在一阶导数过小或初值离根太远的情况下。

因此,使用时需要谨慎,并进行必要的检查和处理。

功能强大的牛顿插值matlab程序代码

功能强大的牛顿插值matlab程序代码

format long;way_in = input('请选择输入的内容(1或2):\n1、输入为f(x)表达式,区间[a,b]及其等分数n的值\n2、输入为f(x)表达式和插值点横坐标xi的值\n');switch way_incase 1f = input('请输入函数表达式:f(x) = ', 's');a = input('请输入区间左端值a:');b = input('请输入区间右端值b:');n = input('请输入区间等分值n:');np = input('请输入插值函数在区间内绘图点数(默认输入100):');for i=1:n+1x(i) = a + (b-a)/n*(i-1);y(i,1) = eval(subs(f,'x(i)','x'));endfor j=1:nfor k=j:ntemp=y(k+1,j)-y(k,j);y(k+1,j+1)=temp/(x(k+1)-x(k+1-j)) ;endc(j)=y(j,j);endc(j+1)=y(j+1,j+1);for k=1:np-1xx(k)= a + (b-a)/np*k;yy(k) = eval(subs(f,'xx(k)','x'));endfor k=1:np-1xs=xx(k);for i=1:n+1if i==1s(i)=c(i);elses(i)=c(i);for j=1:i-1s(i)=s(i)*(xs-x(j));endendendNn(k)=sum(s);way_out = input('请选择要绘出的曲线(1、2或3):\n1、同时输出原始曲线f(x)和插值曲线\n2、只输出插值曲线\n3、只输出原始曲线\n');switch way_outcase 1figure;plot(xx,yy,'r');grid on;hold on;plot(xx,Nn,'b');legend('原始曲线f(x)','插值曲线N(x)');title('牛顿插值');case 2figure;plot(xx,Nn,'m');legend('插值曲线N(x)');title('牛顿插值');case 3figure;plot(xx,yy,'g');legend('原始曲线f(x)');title('牛顿插值');otherwiseerrordlg('请正确选择,输入只能为1、2或者3!','提示','on');endcase 2f = input('请输入函数表达式:f(x) = ', 's');xb = input('请输入插值节点的横坐标x:','s');x = sscanf(xb,'%f');disp('x0,x1,...,xi分别为:');disp(x);n = size(x,1) - 1;if n<2errordlg('请至少输入3个xi的值','提示','on');return;endnp = input('请输入插值函数在区间内绘图点数(默认输入100):');a = x(1);b = x(n+1);for i=1:n+1y(i,1) = eval(subs(f,'x(i)','x'));endfor j=1:nfor k=j:ntemp=y(k+1,j)-y(k,j);y(k+1,j+1)=temp/(x(k+1)-x(k+1-j)) ;endc(j)=y(j,j);endc(j+1)=y(j+1,j+1);for k=1:np-1xx(k) = a + (b-a)/np*k;yy(k) = eval(subs(f,'xx(k)','x'));endfor k=1:np-1xs=xx(k);for i=1:n+1if i==1s(i)=c(i);elses(i)=c(i);for j=1:i-1s(i)=s(i)*(xs-x(j));endendendNn(k)=sum(s);endway_out = input('请选择要绘出的曲线(1、2或3):\n1、同时输出原始曲线f(x)和插值曲线\n2、只输出插值曲线\n3、只输出原始曲线\n');switch way_outcase 1figure;plot(xx,yy,'r');grid on;hold on;plot(xx,Nn,'b');legend('原始曲线f(x)','插值曲线N(x)');title('牛顿插值');case 2figure;plot(xx,Nn,'m');legend('插值曲线N(x)');title('牛顿插值');case 3figure;plot(xx,yy,'g');legend('原始曲线f(x)');title('牛顿插值');otherwiseerrordlg('请正确选择,输入只能为1、2或者3!','提示','on');endotherwiseerrordlg('请正确选择,输入只能为1或2!','提示','on');end。

数值分析matlab代码

数值分析matlab代码

1、%用牛顿法求f(x)=x-sin x 的零点,e=10^(-6)disp('牛顿法');i=1;n0=180;p0=pi/3;tol=10^(-6);for i=1:n0p=p0-(p0-sin(p0))/(1-cos(p0));if abs(p-p0)<=10^(-6)disp('用牛顿法求得方程的根为')disp(p);disp('迭代次数为:')disp(i)break;endp0=p;endif i==n0&&~(abs(p-p0)<=10^(-6))disp(n0)disp('次牛顿迭代后无法求出方程的解')end2、disp('Steffensen加速');p0=pi/3;for i=1:n0p1=0.5*p0+0.5*cos(p0);p2=0.5*p1+0.5*cos(p1);p=p0-((p1-p0).^2)./(p2-2.*p1+p0);if abs(p-p0)<=10^(-6)disp('用Steffensen加速求得方程的根为')disp(p);disp('迭代次数为:')disp(i)break;endp0=p;endif i==n0&&~(abs(p-p0)<=10^(-6))disp(n0)disp('次Steffensen加速后无法求出方程的解')end1、%使用二分法找到方程 600 x^4 -550 x^3 +200 x^2 -20 x -1 =0 在区间[0.1,1]上的根,%误差限为 e=10^-4disp('二分法')a=0.2;b=0.26;tol=0.0001;n0=10;fa=600*(a.^4)-550*(a.^3)+200*(a.^2)-20*a-1;for i=1:n0p=(a+b)/2;fp=600*(p.^4)-550*(p.^3)+200*(p.^2)-20*p-1;if fp==0||(abs((b-a)/2)<tol)disp('用二分法求得方程的根p=')disp(p)disp('二分迭代次数为:')disp(i)break;endif fa*fp>0a=p;else b=p;endendif i==n0&&~(fp==0||(abs((b-a)/2)<tol))disp(n0)disp('次二分迭代后没有求出方程的根')end2、%使用牛顿法找到方程 600 x^4 -550 x^3 +200 x^2 -20 x -1 =0 在区间[0.1,1]上的根,%误差限为 e=10^-4disp('牛顿法')p0=0.3;for i=1:n0p=p0-(600*(p0.^4)-550*(p0.^3)+200*(p0.^2)-20*p0-1)./(2400*(p0.^3) -1650*p0.^2+400*p0-20);if(abs(p-p0)<tol)disp('用牛顿法求得方程的根p=')disp(p)disp('牛顿迭代次数为:')disp(i)break;endp0=p;endif i==n0&&~(abs(p-p0)<tol)disp(n0)disp('次牛顿迭代后没有求出方程的根')end3、%使用割线法找到方程 600 x^4 -550 x^3 +200 x^2 -20 x -1 =0 在区间[0.1,1]上的根,%误差限为 e=10^-4disp('割线法')p0=0.2;p1=0.25;q0=600*(p0.^4)-550*(p0.^3)+200*(p0.^2)-20*p0-1;q1=600*(p1.^4)-550*(p1.^3)+200*(p1.^2)-20*p1-1;for i=2:n0p=p1-q1*(p1-p0)/(q1-q0);if abs(p-p1)<toldisp('用割线法求得方程的根p=')disp(p)disp('割线法迭代次数为:')disp(i)break;endp0=p1;q0=q1;pp=p1;p1=p;q1=600*(p.^4)-550*(p.^3)+200*(p.^2)-20*p-1;endif i==n0&&~(abs(p-pp)<tol)disp(n0)disp('次割线法迭代后没有求出方程的根')end4、%使用试位法找到方程 600 x^4 -550 x^3 +200 x^2 -20 x -1 =0 在区间[0.1,1]上的根,%误差限为 e=10^-4disp('试位法')p0=0.2;p1=0.25;q0=600*(p0.^4)-550*(p0.^3)+200*(p0.^2)-20*p0-1;q1=600*(p1.^4)-550*(p1.^3)+200*(p1.^2)-20*p1-1;for i=2:n0p=p1-q1*(p1-p0)/(q1-q0);if abs(p-p1)<toldisp('用试位法求得方程的根p=')disp(p)disp('试位法迭代次数为:')disp(i)break;endq=600*(p.^4)-550*(p.^3)+200*(p.^2)-20*p-1;if q*q1<0p0=p1;q0=q1;endpp=p1;p1=p;q1=q;endif i==n0&&~(abs(p-pp)<tol)disp(n0)disp('次试位法迭代后没有求出方程的根')end5、%使用muller方法找到方程 600 x^4 -550 x^3 +200 x^2 -20 x -1 =0 在区间[0.1,1]上的根,%误差限为 e=10^-4disp('muller法')x0=0.1;x1=0.2;x2=0.25;h1=x1-x0;h2=x2-x1;d1=((600*(x1.^4)-550*(x1.^3)+200*(x1.^2)-20*x1-1)-(600*(x0.^4)-55 0*(x0.^3)+200*(x0.^2)-20*x0-1))/h1;d2=((600*(x2.^4)-550*(x2.^3)+200*(x2.^2)-20*x2-1)-(600*(x1.^4)-55 0*(x1.^3)+200*(x1.^2)-20*x1-1))/h2;d=(d2-d1)/(h2+h1);for i=3:n0b=d2+h2*d;D=(b*b-4*(600*(x2.^4)-550*(x2.^3)+200*(x2.^2)-20*x2-1)*d)^0.5;if(abs(d-D)<abs(d+D))E=b+D;else E=b-D;endh=-2*(600*(x2.^4)-550*(x2.^3)+200*(x2.^2)-20*x2-1)/E;p=x2+h;if abs(h)<toldisp('用muller方法求得方程的根p=')disp(p)disp('muller方法迭代次数为:')disp(i)break;endx0=x1;x1=x2;x2=p;h1=x1-x0;h2=x2-x1;d1=((600*(x1.^4)-550*(x1.^3)+200*(x1.^2)-20*x1-1)-(600*(x0.^4)-55 0*(x0.^3)+200*(x0.^2)-20*x0-1))/h1;d2=((600*(x2.^4)-550*(x2.^3)+200*(x2.^2)-20*x2-1)-(600*(x1.^4)-55 0*(x1.^3)+200*(x1.^2)-20*x1-1))/h2;d=(d2-d1)/(h2+h1);endif i==n0%条件有待商榷?!disp(n0)disp('次muller方法迭代后没有求出方程的根')end1、%观察Lagrange插值的Runge现象x=-1:0.05:1;y=1./(1+25.*x.*x);plot(x,y),grid on;n=5;x=-1:2/n:1;y=1./(1+25.*x.*x);for i=1:n+1q(1,i)=y(i);endh=0.05;z=-1:h:1;for k=1:2/h+1for i=2:n+1for j=2:iq(j,i)=((z(k)-x(i-j+1))*q(j-1,i)-(z(k)-x(i))*q(j-1,i-1))/(x(i)-x( i-j+1));endendw(k)=q(n+1,n+1);endhold on, plot(z,w,'r'),grid on;%**** n=10 ****n=10;x=-1:2/n:1;y=1./(1+25.*x.*x);for i=1:n+1q(1,i)=y(i);endh=0.05;z=-1:h:1;for k=1:2/h+1for i=2:n+1for j=2:iq(j,i)=((z(k)-x(i-j+1))*q(j-1,i)-(z(k)-x(i))*q(j-1,i-1))/(x(i)-x( i-j+1));endendw(k)=q(n+1,n+1);endhold on,plot(z,w,'k'),grid on;legend ('原始图','n=5','n=10');2、%固支样条插植%********第一段********x=[1,2,5,6,7,8,10,13,17];a=[3,3.7,3.9,4.2,5.7,6.6,7.1,6.7,4.5];n=numel(a);for i=1:n-1h(i)=x(i+1)-x(i);endA=[2*h(1),h(1),0,0,0,0,0,0,0;h(1),2*(h(1)+h(2)),h(2),0,0,0,0,0,0;0,h(2),2*(h(2)+h(3)),h(3),0,0,0,0,0;0,0,h(3),2*(h(3)+h(4)),h(4),0,0,0,0;0,0,0,h(4),2*(h(4)+h(5)),h(5),0,0,0;0,0,0,0,h(5),2*(h(5)+h(6)),h(6),0,0;0,0,0,0,0,h(6),2*(h(6)+h(7)),h(7),0;0,0,0,0,0,0,h(7),2*(h(7)+h(8)),h(8);0,0,0,0,0,0,0,h(8),2*h(8)];e=[3*(a(2)-a(1))/h(1)-3;3*(a(3)-a(2))/h(2)-3*(a(2)-a(1))/h(1);3*(a(4)-a(3))/h(3)-3*(a(3)-a(2))/h(2);3*(a(5)-a(4))/h(4)-3*(a(4)-a(3))/h(3);3*(a(6)-a(5))/h(5)-3*(a(5)-a(4))/h(4);3*(a(7)-a(6))/h(6)-3*(a(6)-a(5))/h(5);3*(a(8)-a(7))/h(7)-3*(a(7)-a(6))/h(6);3*(a(9)-a(8))/h(8)-3*(a(8)-a(7))/h(7);3*(-0.67)-3*(a(9)-a(8))/h(8)];c=inv(A)*e;for i=1:8b(i)=(a(i+1)-a(i))/h(i)-h(i)*(2*c(i)+c(i+1))/3;d(i)=(c(i+1)-c(i))/(3*h(i));endfor i=1:8z=x(i):0.05:x(i+1);w=a(i)+b(i).*(z-x(i))+c(i).*(z-x(i)).^2+d(i).*(z-x(i)).^3; grid on, plot(z,w),hold on;end%********第二段********x=[17,20,23,24,25,27,27.7];a=[4.5,7,6.1,5.6,5.8,5.2,4.1];for i=1:6h(i)=x(i+1)-x(i);endA=[2*h(1),h(1),0,0,0,0,0;h(1),2*(h(1)+h(2)),h(2),0,0,0,0;0,h(2),2*(h(2)+h(3)),h(3),0,0,0;0,0,h(3),2*(h(3)+h(4)),h(4),0,0;0,0,0,h(4),2*(h(4)+h(5)),h(5),0;0,0,0,0,h(5),2*(h(5)+h(6)),h(6)0,0,0,0,0,h(6),2*h(6)];e=[3*(a(2)-a(1))/h(1)-3*3;3*(a(3)-a(2))/h(2)-3*(a(2)-a(1))/h(1);3*(a(4)-a(3))/h(3)-3*(a(3)-a(2))/h(2);3*(a(5)-a(4))/h(4)-3*(a(4)-a(3))/h(3);3*(a(6)-a(5))/h(5)-3*(a(5)-a(4))/h(4);3*(a(7)-a(6))/h(6)-3*(a(6)-a(5))/h(5);3*(-4)-3*(a(7)-a(6))/h(6)];c=inv(A)*e;for i=1:6b(i)=(a(i+1)-a(i))/h(i)-h(i)*(2*c(i)+c(i+1))/3;d(i)=(c(i+1)-c(i))/(3*h(i));endfor i=1:6z=x(i):0.05:x(i+1);w=a(i)+b(i).*(z-x(i))+c(i).*(z-x(i)).^2+d(i).*(z-x(i)).^3; grid on, plot(z,w),hold on;end%********第三段********x=[27.7,28,29,30];a=[4.1,4.3,4.1,3];for i=1:3h(i)=x(i+1)-x(i);endA=[2*h(1),h(1),0,0;h(1),2*(h(1)+h(2)),h(2),0;0,h(2),2*(h(2)+h(3)),h(3);0,0,h(3),2*h(3)];e=[3*(a(2)-a(1))/h(1)-3*0.33;3*(a(3)-a(2))/h(2)-3*(a(2)-a(1))/h(1);3*(a(4)-a(3))/h(3)-3*(a(3)-a(2))/h(2);3*(-1.5)-3*(a(4)-a(3))/h(3)];c=inv(A)*e;for i=1:3b(i)=(a(i+1)-a(i))/h(i)-h(i)*(2*c(i)+c(i+1))/3;d(i)=(c(i+1)-c(i))/(3*h(i));endfor i=1:3z=x(i):0.05:x(i+1);w=a(i)+b(i).*(z-x(i))+c(i).*(z-x(i)).^2+d(i).*(z-x(i)).^3; grid on, plot(z,w),hold on;endgrid on,title('注:横纵坐标的比例不一样!!!');1、%用不动点迭代法求方程 x-e^x+4=0的正根与负根,误差限是10^-6%disp('不动点迭代法');n0=100;p0=-5;for i=1:n0p=exp(p0)-4;if abs(p-p0)<=10^(-6)if p<0disp('|p-p0|=')disp(abs(p-p0))disp('不动点迭代法求得方程的负根为:')disp(p);break;elsedisp('不动点迭代法无法求出方程的负根.')endelsep0=p;endendif i==n0disp(n0)disp('次不动点迭代后无法求出方程的负根')endp1=1.7;for i=1:n0pp=exp(p1)-4;if abs(pp-p1)<=10^(-6)if pp>0disp('|p-p1|=')disp(abs(pp-p1))disp('用不动点迭代法求得方程的正根为')disp(pp);elsedisp('用不动点迭代法无法求出方程的正根');endbreak;elsep1=pp;endendif i==n0disp(n0)disp('次不动点迭代后无法求出方程的正根')end2、%用牛顿法求方程 x-e^x+4=0的正根与负根,误差限是10^-6 disp('牛顿法')n0=80;p0=1;for i=1:n0p=p0-(p0-exp(p0)+4)/(1-exp(p0));if abs(p-p0)<=10^(-6)disp('|p-p0|=')disp(abs(p-p0))disp('用牛顿法求得方程的正根为')disp(p);break;elsep0=p;endendif i==n0disp(n0)disp('次牛顿迭代后无法求出方程的解')endp1=-3;for i=1:n0p=p1-(p1-exp(p1)+4)/(1-exp(p1));if abs(p-p1)<=10^(-6)disp('|p-p1|=')disp(abs(p-p1))disp('用牛顿法求得方程的负根为')disp(p);break;elsep1=p;endendif i==n0disp(n0)disp('次牛顿迭代后无法求出方程的解')end1、使用欧拉法、改进欧拉法和四阶R-K方法求下列微分方程的解。

插值法

插值法

实验四 插值法一、题目:将12)(26++=x x x f 在【0,2】上按5个等距插值基点确定的数值作为插值原则,给出Newton 插值多项式,计算)3.1(n N ,并计算余项)3.1(n R 。

用matlab 画出两条曲线的图形二、方法:Newton 插值法三、程序:clearclcx0=0:0.4:2;y0=x0.^6+2.*x0.^2+1;syms xfor i=1:5a=1;for j=1:5if j~=ia=expand(a*(x-x0(j)));endendb=1;for k=1:5if k~=ib=b*(x0(i)-x0(k));endendA(i)=expand(a/b);endN=0;for p=1:5N=N+y0(p)*A(p);endN四、结果:N =5620492334958391409/540431955284459520*x^4-10808639105689222577/562949953421312000*x^3+239573485777601629553/16888498602639360000*x^2-17293822569102771343/7036874417766400000*x+1>> x=1.3;>>5620492334958391409/540431955284459520*x^4-10808639105689222577/562949953421312000*x^3+239573485777601629553/16888498602639360000*x^2-17293822569102771343/7036874417766400000*x+1ans =9.29982400000001>>x.^6+2.*x.^2+1-(5620492334958391409/540431955284459520*x^4-10808639105689222577/56 2949953421312000*x^3+239573485777601629553/16888498602639360000*x^2-17293822569 102771343/7036874417766400000*x+1)ans =-0.09301500000001N9.29982400000001)3.1(nR(1.3)= -0.09301500000001n画图:>> x1=0:0.0001:1;>> y1=x1.^6+2.*x1.^2+1;>>y2=5620492334958391409/540431955284459520*x1.^4-10808639105689222577/56294995342 1312000*x1.^3+239573485777601629553/16888498602639360000*x1.^2-17293822569102771 343/7036874417766400000*x1+1;>> plot(x1,y1,'r-')>> hold on>> plot(x1,y2,'k-')五、拓展:牛顿插值法有很好的承袭性,增加或减少节点时只增加或减少一行数据,有利于计算。

matlab编程实现二分法,牛顿法,黄金分割法,最速下降matlab程序代码

matlab编程实现二分法,牛顿法,黄金分割法,最速下降matlab程序代码

matlab编程实现二分法,牛顿法,黄金分割法,最速下降matlab程序代码用二4224min ()f t t t t =--[,.]t ∈内的极小值点,要求准1.function [t d]=erfenfa(a,b)k=1; %记录循环次数 while abs(a-b)>0.0005c=(a+b)/2;C(k)=c; %存储每次循环中点c 的值if ff(c)<0a=c;endif ff(c)==0t1=c;break ;endif ff(c)>0b=c;endk=k+1;endt=(a+b)/2; %最终符合要求的值d=f(t); %最优解Ckfunction y=f(t)y=t^4-2*t^2-4*t;function y=ff(t)y=4*t^3-4*t-4;运行结果>> [t d]=erfenfa(1,1.5)C =Columns 1 through 91.2500 1.3750 1.3125 1.3438 1.3281 1.3203 1.3242 1.3262 1.3252Column 101.3247k =11t =1.3250d =-5.72902.黄金分割法 f (x)=x3-2x+1 初始区间[0, 3],收敛精度0.5 function [t,f]=huangjinfenge(a,b)m=1-(sqrt(5)-1)/2;t2=a+m*(b-a)f2=g(t2);t1=a+b-t2f1=g(t1);while abs(t1-t2)>0.5if f1<f2< bdsfid="121" p=""></f2<>a=t2;t2=t1f2=f1;t1=a+b-t2f1=g(t1);elseb=t1;t1=t2f1=f2;t2=a+m*(b-a)f2=g(t2);endendt=(t1+t2)/2;f=g(t);function y=g(t)y=t^3-2*t+1;运行结果> [t,f]=huangjinfenge(0,3)t2 =1.1459t1 =1.8541t1 =1.1459t2 =0.7082t =0.9271f =-0.0574>>3. 用牛顿法求解291min ()sin f x x x =--初始迭代点为x 0=0.4, 要求准确到小数点后第5位小数function [t1,d]=Newton(t0)t=t0-ff(t0)/fff(t0);k=1;%记录迭代次数T(1)=t;%存储迭代点while abs(t-t0)>0.000005t0=t;t=t0-ff(t)/fff(t);k=k+1;T(k)=t;endt1=t0;d=f(t1);kTfunction y=f(x)y=9*x^2-sin(x)-1;function y=ff(x)y=18*x-cos(x);function y=fff(x)y=18+sin(x);运行结果>> [t1,d]=Newton(0.4)k =3T =0.0586 0.0555 0.0555t1 =0.0555d =-1.0277>>4. 最速下降法验证课本上的例题求解291min ()sin f x x x =--初始迭代点为x 0=0.4, 要求准确到小数点后第5位小数function [G,g,X,F]=zuisu(X0)F(1)=f(X0);%存储x 点处的值G(:,1)=h(X0); %存储梯度向量g(1)=norm(G(:,1));%存储梯度模长X(:,1)=X0; %存储x 值A=[2,0;0,8];for j=1:2X(:,j+1)=X(:,j)-(G(:,j)'*G(:,j))/(G(:,j)'*A*G(:,j))*G(:,j);F(j+1)=f(X(:,j+1));G(:,j+1)=h(X(:,j+1));g(j+1)=norm(G(:,j+1));endif (G(:,2)'*G(:,1)<1E-10& G(:,3)'*G(:,2)<1E-10)disp(['相邻两搜索方向是正交的'])endfunction y=f(X)y=X(1)^2+4*X(2)^2;function n=h(X)n=[2*X(1),8*X(2)]';运行结果>> [G,g,X,F]=zuisu(X0)相邻两搜索方向是正交的G =2.0000 1.4769 0.2215 8.0000 -0.3692 0.8862g =8.2462 1.5224 0.9134X =1.0000 0.7385 0.1108 1.0000 -0.0462 0.1108F =5.0000 0.5538 0.0613 >>。

matlab牛顿插值法函数

matlab牛顿插值法函数

matlab牛顿插值法函数
牛顿插值法是一种常用的数值计算方法,可以通过已知的离散数据点来估计未知的函数值。

该方法基于插值多项式的思想,利用已知数据点的信息来构建一个多项式,然后利用该多项式来估计其他点的函数值。

牛顿插值法的基本思想是利用差商的概念,通过递推的方式来计算插值多项式的系数。

具体来说,给定n个数据点(x0, y0), (x1, y1), ..., (xn, yn),我们可以得到一个n次插值多项式。

该多项式的形式为:
P(x) = f[x0] + (x - x0)f[x0, x1] + (x - x0)(x - x1)f[x0, x1, x2] + ... + (x - x0)(x - x1)...(x - xn-1)f[x0, x1, ..., xn]
其中,f[xi]表示差商,f[x0, x1]表示二阶差商,以此类推。

牛顿插值法的优点是计算简单,且对于多项式插值问题具有很高的精度。

然而,该方法也有一些局限性。

首先,插值多项式的次数随着数据点的增加而增加,可能导致多项式振荡或者过拟合。

此外,当数据点不均匀分布时,插值多项式的精度可能会受到影响。

为了解决这些问题,可以使用其他插值方法,如拉格朗日插值法或样条插值法。

这些方法在一定程度上克服了牛顿插值法的局限性,但也引入了一些新的问题。

因此,在实际应用中,需要根据具体问题的特点选择合适的插值方法。

牛顿插值法是一种常用的数值计算方法,可以用于估计未知函数值。

虽然该方法有一些局限性,但在适当的条件下,仍然可以获得较高的插值精度。

在实际应用中,需要综合考虑问题的特点和需求,选择合适的插值方法。

拉格朗日差值、牛顿插值以及三次样条插值的matlab实现

拉格朗日差值、牛顿插值以及三次样条插值的matlab实现

拉格朗日差值、牛顿插值以及三次样条插值的matlab实现% Lagrange插值clearclc%-----------------------------n=10; %结点个数lb=-1; %下界ub=1; %上界step=0.01; %作图点步长%-----------------------------% 原始函数图形x0=lb:step:ub;y0=1./(1+25*x0.^2);plot(x0,y0,'r-');hold on%-----------------------------% 插值函数for i=1:n+1xi(i)=lb+(ub-lb)*(i-1)/n;yi(i)=1/(1+25*xi(i)^2);end%------------------------------count=1;for x=lb:step:ubfl=0;%--------------------------%求出pn(xk)for k=1:n+1up=1;dn=1;%----------------------%求出f(xk)for i=1:n+1if k~=iup=up*(x-xi(i));dn=dn*(xi(k)-xi(i));endend%----------------------fl=fl+yi(k)*up/dn;endpn(count)=fl;%--------------------------fi(count)=1/(1+25*x^2);%求原函数的值count=count+1;end%------------------------------% L插值函数图x=lb:step:ub;plot(x,pn,'g--')%------------------------------num=(ub-lb)/step+1;for i=1:nump_f(i)=pn(i)-fi(i);endcenter=fix(num/2);scale=fix(num/10);a=center-scale;b=center+scale;disp ' pn(i)-fi(i) 的值为:'p_f(a:b) %%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%% %%clear allclc%Newton迭代法求解极小值点%===================================== disp '几点说明:'disp '1.程序中的函数采用课本P102例3.3.2。

  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
相关文档
最新文档