直角坐标牛顿-拉夫逊法潮流计算matlab程序(仅供参考)
%该程序仅针对《电力系统分析下》何仰赞P61 的4节点算例。
%节点电压用直角坐标表示时的牛顿-拉夫逊法潮流计算（matlab程序，可能有小错误，仅供学习之用，如果想要通用的程序，可自己动手改或再到pudn、csdn等网站搜索更好的）%南昌大学电力061 李圣涛2009年5月编写，2012年5月上传
clc %清空command windows
clear all %清空workspace
%为了提高可移植性、可读性、通用性，设置以下变量
N=4; %独立节点数
NPQ=2; %PQ节点数
NPV=1; %PV节点数
K=0; %迭代次数
%请输入最大迭代次数Kmax。可从0开始，以观察第Kmax次迭代的结果
Kmax=input('\n\n 请输入最大迭代次数后回车(可从零开始) Kmax=\n');
small=10^(-5); %ε不能太小
%i为节点标号，其中1号……NPQ号为PQ节点，（NPQ+1）
%号……（N-1）号为PV节点，N号节点为平衡节点
%节点导纳矩阵Y的实部
G=[1.042093 -0.588235 0 -0.453858;
-0.588235 1.069005 0 -0.480769;
0 0 0 0 ;
-0.453858 -0.480769 0 0.934627 ];
%节点导纳矩阵Y的虚部
B=[ -8.242876 2.352941 3.666667 1.891074 ;
2.352941 -4.727377 0 2.403846 ;
3.666667 0 -3.3333333 0 ;
1.891074
2.403846 0 -4.261590 ];
%Y矩阵
Y=complex(G,B);
%给定PQ节点的Pnode、Qnode，PV节点的Pnode、Vnode。（Vnode为节点电压的幅值）Pnode=[ -0.3 -0.55 0.5 ];%PQ、PV节点的初值P
Qnode=[ -0.18 -0.13 0 ];%PQ节点的初值Q
Vnode=[ 0 0 1.10 ];%PV节点的初值V
%迭代初值
e=[ 1.0 1.0 1.1 1.05 ];
f=[ 0 0 0 0 ];
%利用for循环来实现多次迭代。
for K=0:Kmax,
%计算△W。△W=【△P1,△Q1,△P2,△Q2,……,△P,△V^2】(平衡节点不参与迭代) %用公式（11-46）计算NPQ个PQ节点的△P、△Q。算法：先初始化，再不断地累减
%初始化,得△P=【Pnode(1),Pnode(2),……，Pnode(NPQ)】
for i=1:NPQ,
dP(i)=Pnode(i);
dQ(i)=Qnode(i);
end
%不断地累减,得△P和△Q
for i=1:NPQ,
for j=1:N,
dP(i)=dP(i)-e(i)*G(i,j)*e(j)+e(i)*B(i,j)*f(j)-f(i)*G(i,j)*f(i)-f(i)*B(i,j)*e(j);%△P
dQ(i)=dQ(i)-f(i)*G(i,j)*e(j)+f(i)*B(i,j)*f(j)+e(i)*G(i,j)*f(i)+e(i)*B(i,j)*e(j);%△Q
end
end
%用公式（11-47）计算NPV个PV节点的△P、（△V）^2 (其中△P的计算同上)算法：先初始化，再不断地累减
%初始化
for i=(NPQ+1):(N-1),
dP(i)=Pnode(i);
end
%不断地累减,得△P和（△V）^2
for i=(NPQ+1):(N-1),
for j=1:N,
dP(i)=dP(i)-e(i)*G(i,j)*e(j)+e(i)*B(i,j)*f(j)-f(i)*G(i,j)*f(i)-f(i)*B(i,j)*e(j);
dV2(i)=Vnode(i)^2-e(i)^2-f(i)^2;
end
end
%组合成△W=【△P1,△Q1,△P2,△Q2,……,△P,△V^2】(平衡节点不参与迭代)
%将△P间隔地赋入△W
a=1;
for i=1:(N-1),
dW(a)=dP(i);
a=a+2;
end
%将△Q间隔地赋入△W
a=2;
for i=1:NPQ,
dW(a)=dQ(i);
a=a+2;
end
%将△V^2间隔地赋入△W
a=NPQ*2+2;
for i=(NPQ+1):(NPQ+NPV),
dW(a)=dV2(i);
a=a+2;
end
%判断是否小于ε
if max(dW) fprintf('\n 迭代是收敛的。第%d次迭代后终止迭代。\n',K-1); break; end %计算雅克比矩阵各元素。算法：先初始化，再不断地累减 J=zeros(N-1); %初始化 for i=1:N-1, for j=1:N-1, if i<=NPQ %求雅克比矩阵中PQ节点的元素 if i==j J(2*i-1,2*j-1) = -G(i,i)*e(i)-B(i,i)*f(i); J(2*i-1,2*j ) = B(i,i)*e(i)-G(i,i)*f(i); J(2*i ,2*j-1) = B(i,i)*e(i)-G(i,i)*f(i); J(2*i ,2*j ) = G(i,i)*e(i)+B(i,i)*f(i); elseif i~=j J(2*i-1,2*j-1) = -G(i,j)*e(i)-B(i,j)*f(i); J(2*i-1,2*j ) = B(i,j)*e(i)-G(i,j)*f(i); J(2*i ,2*j-1) = B(i,j)*e(i)-G(i,j)*f(i); J(2*i ,2*j ) = G(i,j)*e(i)+B(i,j)*f(i); end end if i>NPQ %求雅克比矩阵中PV节点的元素 if i==j J(2*i-1,2*j-1) = -G(i,i)*e(i)-B(i,i)*f(i); J(2*i-1,2*j ) = B(i,i)*e(i)-G(i,i)*f(i); J(2*i ,2*j-1) = -2*e(i); J(2*i ,2*j ) = -2*f(i); elseif i~=j J(2*i-1,2*j-1) = -G(i,j)*e(i)-B(i,j)*f(i); J(2*i-1,2*j ) = B(i,j)*e(i)-G(i,j)*f(i); J(2*i ,2*j-1) = 0; J(2*i ,2*j ) = 0; end end end end %不断地累减 for i=1:NPQ, for k=1:N, J(2*i-1,2*i-1) = J(2*i-1,2*i-1)-G(i,k)*e(k)+B(i,k)*f(k); J(2*i-1,2*i ) = J(2*i-1,2*i )-G(i,k)*f(k)-B(i,k)*e(k); J(2*i ,2*i-1) = J(2*i ,2*i-1)+G(i,k)*f(k)+B(i,k)*e(k); J(2*i ,2*i ) = J(2*i ,2*i )-G(i,k)*e(k)+B(i,k)*f(k); end end %根据△W=-J*△V, 得△V dV=(-J)\dW'; %用dV对e、f进行修正,并得到复数表示的V for i=1:(N-1), e(i)=e(i)+dV(2*i-1); f(i)=f(i)+dV(2*i ); V(i)=complex(e(i),f(i)); end V(N)=complex(e(N),f(N)); format long g; fprintf('\n 第%d次迭代(共%d个独立节点)',K,N); V end 牛顿插值法 插值法是利用函数f (x)在某区间中若干点的函数值，作出适当的特定函数，在这些点上取已知值，在区间的其他点上用这特定函数的值作为函数f (x)的近似值。如果这特定函数是多项式，就称它为插值多项式。当插值节点增减时全部插值基函数均要随之变化，这在实际计算中很不方便。为了克服这一缺点，提出了牛顿插值。牛顿插值通过求各阶差商，递推得到的一个公式： f(x)=f[x0]+f[x0,x1](x-x0)+f[x0,x1,x2](x-x0)(x-x1)+...f[x0,...xn](x-x0 )...(x-xn-1)+Rn(x)。 插值函数 插值函数的概念及相关性质[1] 定义：设连续函数y-f(x) 在区间[a,b]上有定义，已知在n+1个互异的点 x0,x1,…xn上取值分别为y0,y1,…yn （设a≤ x1≤x2……≤xn≤b)。若在函数类中存在以简单函数P(x) ，使得P(xi)=yi,则称P(x) 为f(x)的插值函数. 称x1,x2,…xn 为插值节点，称[a,b]为插值区间。 定理：n次代数插值问题的解存在且唯一。 牛顿插值法C程序 程序框图#include 牛顿插值法在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页的表 d(:,1)=y ’; %此句是把向量y 的转置，即123()()()f x y f x f x ?? ?= ? ?? ?，赋值给零矩阵d 的第一列 %下面运用两个for 循环来构造书中21页的差商表 %第一个循环（父循环），循环变量为k for k=2:n %用来表示零矩阵d 中的第几行 %第二个循环（父循环），循环变量为k for 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 -，这种计算差商的方法是根据差商表的排列位置而得来，具体解释见下面。 end end %下面以二次牛顿插值公式为例解析双循环构造差商表，让我们先来看看构造好的差商表 121232312333 () () [,] ()[,][,,]X f x d f x f x x f x f x x f x x x ????=?????? 前言 电力系统潮流计算是研究电力系统稳态运行情况的一种计算,它根据给定的运行件及系统接线情况确定整个电力系统各部分的运行状态。在电力系统规划设计和现有电力系统运行方式的研究中,都需要利用潮流计算来定量分析、比较供电方案或运行方式的合理性、可靠性和经济性。本次课程设计任务是闭式网络的潮流计算，用到的方法为牛顿拉夫逊极坐标法潮流计算。 牛顿法是数学中解决非线性方程式的典型方法，有较好的收敛性。解决电力系统潮流计算问题是以导纳距阵为基础的，因此，只要在迭代过程中尽可能保持方程式系数距阵的稀疏性，就可以大大提高牛顿法潮流程序的放率。自从20 世纪60 年代中期利用了最佳顺序消去法以后，牛顿法在收敛性、内存要求、速度方面都超过了阻抗法，成为直到目前仍在广泛采用的优秀方法。 目录 1任务书 (2) 2.模型简介及等值电路 (3) 3.设计原理 (4) 4.修正方程的建立 (7) 5.程序流程图及MATLAB程序编写 (9) 6.结果分析 (15) 7.设计总结 (25) 8.参考文献 (26) 《电力系统分析》 课程设计任务书 题目极坐标表示的牛顿拉夫逊法潮流计算程序设计学生姓名学号专业班级 设计内容与要求1. 设计要求 掌握MATLAB语言编程方法；理解和掌握运用计算机进行潮流计算的基本算法原理；针对某一具体电网，进行潮流计算程序设计。 其目的在于加深学生对电力系统稳态分析中课程中基本概念和计算方法的理解，培养学生运用所学知识分析和解决问题的能力。 2. 内容 1）学习并掌握MATLAB语言。 2）掌握变压器非标准变比概念及非标准变比变压器的等值电路。掌握节点导纳矩阵的概念及导纳矩阵的形成和修改方法。 3）掌握电力系统功率方程、变量和节点分类。 4）掌握利用极坐标表示的牛-拉法进行潮流计算的方法和步骤。 5）选择一个某一具体电网，编制程序流程框图。 6）利用MATLAB语言编写该模型的潮流计算程序，并上机调试程序，对计算结果进行分析。 7）整理课程设计论文。 起止时间2013 年7 月 4 日至2013 年7月10日指导教师签名年月日 系（教研室）主任签 名 年月日学生签名年月日 2 模型简介及等值电路 2.1课程设计模型：模型3 Matlab实现数值分析插值及积分 摘要： 数值分析(numerical analysis)是研究分析用计算机求解数学计算问题的数值计算方法及其理论的学科，是数学的一个分支，它以数字计算机求解数学问题的理论和方法为研究对象。在实际生产实践中，常常将实际问题转化为数学模型来解决，这个过程就是数学建模。学习数值分析这门课程可以让我们学到很多的数学建模方法。 分别运用matlab数学软件编程来解决插值问题和数值积分问题。题目中的要求是计算差值和积分，对于问题一，可以分别利用朗格朗日插值公式，牛顿插值公式，埃特金逐次线性插值公式来进行编程求解，具体matlab代码见正文。编程求解出来的结果为：=+。 其中Aitken插值计算的结果图如下： 对于问题二，可以分别利用复化梯形公式，复化的辛卜生公式，复化的柯特斯公式编写程序来进行求解，具体matlab代码见正文。编程求解出来的结果为： 0.6932 其中复化梯形公式计算的结果图如下： 问题重述 问题一：已知列表函数 表格 1 分别用拉格朗日，牛顿，埃特金插值方法计算。 问题二：用复化的梯形公式，复化的辛卜生公式，复化的柯特斯公式计算积分，使精度小于5。 问题解决 问题一：插值方法 对于问题一，用三种差值方法：拉格朗日，牛顿，埃特金差值方法来解决。 一、拉格朗日插值法： 拉格朗日插值多项式如下： 首先构造1+n 个插值节点n x x x ,,,10 上的n 插值基函数，对任一点i x 所对应的插值基函数 )(x l i ，由于在所有),,1,1,,1,0(n i i j x j +-=取零值，因此)(x l i 有因子 )())(()(110n i i x x x x x x x x ----+- 。又因)(x l i 是一个次数不超过n 的多项式，所以只 可能相差一个常数因子，固)(x l i 可表示成： )())(()()(110n i i i x x x x x x x x A x l ----=+- 利用1)(=i i x l 得： 摘要 本文，首先简单介绍了基于在MALAB中行潮流计算的原理、意义，然后用具体的实例，简单介绍了如何利用MALAB去进行电力系统中的潮流计算。 众所周知，电力系统潮流计算是研究电力系统稳态运行情况的一种计算，它根据给定的运行条件及系统接线情况确定整个电力系统各部分的运行状态：各线的电压、各元件中流过的功率、系统的功率损耗等等。在电力系统规划的设计和现有电力系统运行方式的研究中，都需要利用潮流计算来定量地分析比较供电方案或运行方式的合理性、可靠性和经济性。 此外，在进行电力系统静态及暂态稳定计算时，要利用潮流计算的结果作为其计算的基础；一些故障分析以及优化计算也需要有相应的潮流计算作配合；潮流计算往往成为上述计算程序的一个重要组成部分。以上这些，主要是在系统规划设计及运行方式安排中的应用，属于离线计算范畴。 牛顿－拉夫逊法在电力系统潮流计算的常用算法之一，它收敛性好，迭代次数少。本文介绍了电力系统潮流计算机辅助分析的基本知识及潮流计算牛顿－拉夫逊法，最后介绍了利用MTALAB程序运行的结果。 关键词：电力系统潮流计算，牛顿－拉夫逊法，MATLAB ABSTRACT This article first introduces the flow calculation based on the principle of MALAB Bank of China, meaning, and then use specific examples, a brief introduction, how to use MALAB to the flow calculation in power systems. As we all know, is the study of power flow calculation of power system steady-state operation of a calculation, which according to the given operating conditions and system wiring the entire power system to determine the operational status of each part: the bus voltage flowing through the components power, system power loss and so on. In power system planning power system design and operation mode of the current study, are required to quantitatively calculated using the trend analysis and comparison of the program or run mode power supply reasonable, reliability and economy. In addition, during the power system static and transient stability calculation, the results of calculation to take advantage of the trend as its basis of calculation; number of fault analysis and optimization also requires a corresponding flow calculation for cooperation; power flow calculation program often become the an important part. These, mainly in the way of system design and operation arrangements in the application areas are off-line calculation. Newton - Raphson power flow calculation in power system is one commonly used method, it is good convergence of the iteration number of small, introduce the trend of computer-aided power system analysis of the basic knowledge and power flow Newton - Raphson method, introduced by the last matlab run results. Keywords：power system flow calculation, Newton – Raphson method, matlab xxxx实验报告 学生姓名：学号：专业班级： 实验类型：□验证□综合■设计□创新实验日期：实验成绩： 一、实验项目名称 P-Q分解法潮流计算实验 二、实验目的与要求： 目的：电力系统分析的潮流计算是电力系统分析的一个重要的部分。通过对电力系统潮流分布的分析和计算，可进一步对系统运行的安全性，经济性进行分析、评估，提出改进措施。电力系统潮流的计算和分析是电力系统运行和规划工作的基础。 潮流计算是指对电力系统正常运行状况的分析和计算。通常需要已知系统参数和条件，给定一些初始条件，从而计算出系统运行的电压和功率等；潮流计算方法很多：高斯-塞德尔法、牛顿-拉夫逊法、P-Q分解法、直流潮流法，以及由高斯-塞德尔法、牛顿-拉夫逊法演变的各种潮流计算方法。 本实验采用P-Q分解法进行电力系统分析的潮流计算程序的编制与调试，获得电力系统中各节点电压，为进一步进行电力系统分析作准备。通过实验教学加深学生对电力系统潮流计算原理的理解和计算，初步学会运用计算机知识解决电力系统的问题，掌握潮流计算的过程及其特点。熟悉各种常用应用软件，熟悉硬件设备的使用方法，加强编制调试计算机程序的能力，提高工程计算的能力，学习如何将理论知识和实际工程问题结合起来。 要求：编制调试电力系统潮流计算的计算机程序。程序要求根据已知的电力网的数学模型（节点导纳矩阵）及各节点参数，完成该电力系统的潮流计算，要求计算出节点电压、功率等参数。 三、主要仪器设备及耗材 每组计算机1台、相关计算软件1套 四、实验内容： 1.理论分析： P-Q分解法是从改进和简化牛顿法潮流程序的基础上提出来的，它的基本思想是：把节点功率表示为电压向量的极坐标方程式，抓住主要矛盾，以有功功率误差作为修正电压向量角度的依据，以无功功率误差作为修正电压幅值的依据，把有功功率和无功功率迭代分开来进行。 牛顿法潮流程序的核心是求解修正方程式，当节点功率方程式采取极坐标系统时， 实验五多项式插值逼近 信息与计算科学金融崔振威201002034031 一、实验目的： 拉格朗日插值和牛顿插值的数值实现 二、实验内容：p171.1、p178.1、龙格现象数值实现 三、实验要求： 1、根据所给题目构造相应的插值多项式， 2、编程实现两类插值多项式的计算 3、试分析多项式插值造成龙格现象的原因 主程序 1、拉格朗日 function [c,l]=lagran(x,y) %c为多项式函数输出的系数 %l为矩阵的系数多项式 %x为横坐标上的坐标向量 %y为纵坐标上的坐标向量 w=length(x); n=w-1; l=zeros(w,w); for k=1:n+1 v=1; for j=1:n+1 if k~=j v=conv(v,poly(x(j)))/(x(k)-x(j)) %对多项式做卷积运算end end l(k,:)=v; end c=y*l; 牛顿插值多项式主程序 function [p2,z]=newTon(x,y,t) %输入参数中x,y为元素个数相等的向量 %t为插入的定点 %p2为所求得的牛顿插值多项式 %z为利用多项式所得的t的函数值。 n=length(x); chaS(1)=y(1); for i=2:n x1=x;y1=y; x1(i+1:n)=[]; y1(i+1:n)=[]; n1=length(x1); s1=0; for j=1:n1 t1=1; for k=1:n1 if k==j %如果相等则跳出循环 continue; else t1=t1*(x1(j)-x1(k)); end end s1=s1+y1(j)/t1; end chaS(i)=s1; end b(1,:)=[zeros(1,n-1) chaS(1)]; cl=cell(1,n-1); %cell定义了一个矩阵 for i=2:n u1=1; for j=1:i-1 u1=conv(u1,[1 -x(j)]); %conv()用于多项式乘法、矩阵乘法 cl{i-1}=u1; end cl{i-1}=chaS(i)*cl{i-1}; b(i,:)=[zeros(1,n-i),cl{i-1}]; end p2=b(1,:); for j=2:n p2=p2+b(j,:); end if length(t)==1 rm=0; for i=1:n rm=rm+p2(i)*t^(n-i); end z=rm; else k1=length(t); rm=zeros(1,k1); for j=1:k1 for i=1:n rm(j)=rm(j)+p2(i)*t(j)^(n-i); end %本程序的功能是用牛顿——拉夫逊法进行潮流计算 % B1矩阵：1、支路首端号；2、末端号；3、支路阻抗；4、支路对地电纳 % 5、支路的变比；6、支路首端处于K侧为1，1侧为0 % B2矩阵：1、该节点发电机功率；2、该节点负荷功率；3、节点电压初始值% 4、PV节点电压V的给定值；5、节点所接的无功补偿设备的容量% 6、节点分类标号：1为平衡节点（应为1号节点）；2为PQ节点；% 3为PV节点； clear; n=10;%input('请输入节点数：n='); nl=10;%input('请输入支路数：nl='); isb=1;%input('请输入平衡母线节点号：isb='); pr=0.00001;%input('请输入误差精度：pr='); B1=[1 2 0.03512+0.08306i 0.13455i 1 0; 2 3 0.0068+0.18375i 0 1.02381 1; 1 4 0.05620+0.13289i 0.05382i 1 0; 4 5 0.00811+0.24549i 0 1.02381 1; 1 6 0.05620+0.13289i 0.05382i 1 0; 4 6 0.04215+0.09967i 0.04037i 1 0; 6 7 0.0068+0.18375i 0 1.02381 1; 6 8 0.02810+0.06645i 0.10764i 1 0; 8 10 0.00811+0.24549i 0 1 1; 8 9 0.03512+0.08306i 0.13455i 1 0] B2=[0 0 1.1 1.1 0 1; 0 0 1 0 0 2; 0 0.343+0.21256i 1 0 0 2; 0 0 1 0 0 2; 0 0.204+0.12638i 1 0 0 2; 0 0 1 0 0 2; 0 0.306+0.18962i 1 0 0 2; 0 0 1 0 0 2; 0.5 0 1.1 1.1 0 3; 0 0.343+0.21256i 1 0 0 2] ;%input('请输入各节点参数形成的矩阵：B2='); Y=zeros(n);e=zeros(1,n);f=zeros(1,n);V=zeros(1,n);sida=zeros(1,n);S1=zeros(nl); % % %--------------------------------------------------- for i=1:nl %支路数 if B1(i,6)==0 %左节点处于1侧 p=B1(i,1);q=B1(i,2); else %左节点处于K侧 p=B1(i,2);q=B1(i,1); end 题目一：多项式插值 某气象观测站在8：00（AM ）开始每隔10分钟对天气作如下观测，用三次多项式插值函数（Newton ）逼近如下曲线，插值节点数据如上表，并求出9点30分该地区的温度（x=10）。 二、数学原理 假设有n+1个不同的节点及函数在节点上的值（x 0，y 0），……（x n ，y n ），插值多项式有如下形式： ）（） ）（（）（）（）（n 10n 102010n x -x )(x -x x -x x P x x x x x x -??-+??+-++=αααα （1） 其中系数i α（i=0,1,2……n ）为特定系数，可由插值样条i i n y x P =） （（i=0,1,2……n ）确定。 根据均差的定义，把x 看成[a,b]上的一点,可得 f(x)= f （0x ）+f[10x x ，]（0x -x ） f[x, 0x ]= f[10x x ，]+f[x,10x x ，] （1x -x ） …… f[x, 0x ，…x 1-n ]= f[x, 0x ，…x n ]+ f[x, 0x ，…x n ](x-x n ) 综合以上式子，把后一式代入前一式，可得到： f(x)= f[0x ]+f[10x x ，]（0x -x ）+ f[210x x x ，，]（0x -x ）（1x -x ）+ …+ f[x, 0x ，…x n ]（0x -x ）…(x-x 1-n )+ f[x, 0x ，…x n ,x ]） （x 1n +ω= N n （x ）+） （x n R 其中 N n （x ）= f[0x ]+f[10x x ，]（0x -x ）+ f[210x x x ，，]（0x -x ）（1x -x ）+ …+ f[x, 0x ，…x n ]（0x -x ）…(x-x 1-n ) （2） ）（x n R = f(x)- N n （x ）= f[x, 0x ， (x) n ,x ]） （x 1n +ω （3） ） （x 1n +ω=（0x -x ）…(x-x n ) Newton 插值的系数i α（i=0,1,2……n ）可以用差商表示。一般有 f k =α[k 10x x x ??，] （k=0,1,2，……，n ） （4） 把（4）代入（1）得到满足插值条件N ）（） （i i n x f x =（i=0,1,2，……n ）的n 次Newton 插值多项式 N n （x ）=f （0x ）+f[10x x ，]（1x -x ）+f[210x x x ，，]（1x -x ）（2x -x ）+……+f[n 10x x x ??，]（1x -x ）（2x -x ）…（1-n x -x ）. 其中插值余项为： ） （）！ （） （）（）（）（x 1n f x N -x f x R 1n 1 n n +++==ωξ ξ介于k 10x x x ??，之间。 三、程序设计 function [y,A,C,L]=newdscg(X,Y,x,M) % y 为对应x 的值，A 为差商表，C 为多项式系数，L 为多项式 % X 为给定节点，Y 为节点值，x 为待求节点 n=length(X); m=length(x); % n 为X 的长度 for t=1:m MATLAB程序设计期中作业 ——编程实现牛顿插值 成员:刘川（P091712797）签名＿＿＿＿＿ 汤意（P091712817）签名＿＿＿＿＿ 王功贺（P091712799）签名＿＿＿＿＿ 班级:2009信息与计算科学 学院:数学与计算机科学学院 日期:2012年05月02日 牛顿插值的算法描述及程序实现 一：问题说明 在我们的实际应用中，通常需要解决这样的问题，通过一些已知的点及其对应的值，去估算另外一些点的值，这些数据之间近似服从一定的规律，于是，这就引入了插值法的思想。 插值法是利用函数f (x)在某区间中若干点的函数值，作出适当的特定函数，在这些点上取已知值，在区间的其他点上用这特定函数的值作为函数f (x)的近似值。如果这特定函数是多项式，就称它为插值多项式。利用插值基函数很容易得到拉格朗日插值多项式，公式结构紧凑，在理论分析中甚为方便，但当插值节点增减时全部插值基函数均要随之变化，整个公式也将发生变化，这在实际计算中是很不方便的，为了克服这一缺点，提出了牛顿插值。 二：算法分析 newton 插值多项式的表达式如下： 010011()()()()()n n n N x c c x x c x x x x x x -=+-+???+--???- 其中每一项的系数c i 的表达式如下： 12011010 [,,,][,,,] [,,,]i i i i i f x x x f x x x c f x x x x x -???-???=???= - 即为f (x)在点01,,,i x x x ???处的i 阶差商，（[]()i i f x f x =，1,2,,i n = ），由差商01[,,,]i f x x x ???的性质可知： () 010 1 [,,,]()i i i j j k j k k j f x x x f x x x ==≠???=-∑∏ 牛顿插值的程序实现方法： 第一步：计算[][][][]001012012,,,,,,,n f x f x x f x x x f x x x x 、、、 、。 第二步：计算牛顿插值多项式中01[,,,]i f x x x ???011()()()i x x x x x x ---???-，1,2,,i n = ，得到n 个多项式。 6.3.5 牛顿插值法的MATLAB 综合程序 求牛顿插值多项式、差商、插值及其误差估计的MATLAB 主程序 function [y,R,A,C,L]=newdscg(X,Y,x,M) n=length(X); m=length(x); for t=1:m z=x(t); A=zeros(n,n);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)));c1=c1*j; 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(k)= polyval(C, z); end R=M*q1/c1;L(k,:)=poly2sym(C); 例6.3.6 给出节点数据00.27)00.4(=-f ，00.1)00.0(=f ，00.2)00.1(=f ，00.17)00.2(=f ，作三阶牛顿插值多项式，计算)345.2(-f ，并估计其误差. 解 首先将名为newdscg.m 的程序保存为M 文件，然后在MATLAB 工作窗口输入程序 >> syms M,X=[-4,0,1,2]; Y =[27,1,2,17]; x=-2.345; [y,R,A,C,P]=newdscg(X,Y,x,M) 运行后输出插值y )345.2(-≈f 及其误差限公式R ，三阶牛顿插值多项式P 及其系数向量C ，差商的矩阵A 如下 y = 22.3211 R = 65133/562949953421312*M （即R =2.3503*M ） A= 27.0000 0 0 0 1.0000 -6.5000 0 0 2.0000 1.0000 1.5000 0 17.0000 15.0000 7.0000 0.9167 C = 0.9167 4.2500 -4.1667 1.0000 P = 11/12*x^3+17/4*x^2-25/6*x+1 (){} 2 1 ()(11),5,10,20: 1252 1()1,(0,1,2,,)()2,(0,1,2,,)() ()2 35,20:1100 (i i i i n n k k k Newton f x x n x f x x i i n f x n x y i n Newton N x S x n x k y f x = -≤≤=+=-+====-+ = 题目：插值多项式和三次样条插值多项式。已知对作、计算函数在点处的值；、求插值数据点 的插值多项式和三次样条插值多项式；、对计算和相应的函数值),()() (1,2,,99)4:()max ()()max ()n k n k n k n k n k n k k k N x S x k E N y N x E S y S x ==-=- 和； 、计算，； 解释你所得到的结果。 算法组织: 本题在算法上需要解决的问题主要是：求出第二问中的Newton 插值多项式 )(x N n 和三次样条插值多项式()n S x 。如此，则第三、四问则迎刃而解。计算两 种插值多项式的算法如下： 一、求Newton 插值多项式)(x N n ，算法组织如下： Newton 插值多项式的表达式如下： )())(()()(110010--???--+???+-+=n n n x x x x x x c x x c c x N 其中每一项的系数c i 的表达式如下： 1102110) ,,,(),,,(),,,(x x x x x f x x x f x x x f c i i i i i -???-???= ???=- 根据i c 以上公式，计算的步骤如下： ?? ??? ?? ?????+??????? ???????????----) ,,,,(1) ,,,(),,,,(),(,),,(2)(,),(),(11101111011010n n n n n n n n x x x x f n x x x f x x x f n x x f x x f x f x f x f 、计算、计算、计算、计算 二、求三次样条插值多项式)(x S n ，算法组织如下： 《工程常用算法》综合实践作业二 完成日期： 2013年 4月 14 日 班级 学号 姓名 主要工作说明 自评成绩 0718 2010071826 崔洪亮 算式与程序的编写 18 0718 2010071815 侯闰上 流程图的编辑，程序的审查 0718 2010071809 赵化川 报告的整理汇总 一.作业题目：三次样条插值与分段插值 已知飞机下轮廓线数据如下： x 3 5 7 9 11 12 13 14 15 y 0 1.2 1.7 2.0 2.1 2.0 1.8 1.2 1.0 1.6 飞机下轮廓线形状大致如下图所示： 要求分别用拉格朗日插值法、Newton 插值法、分段线性插值法和三次样条插值法计算x 每改变0.5时y 的值，即x 取 0.5, 1, 1.5, … , 14.5 时对应的y 值。比较采用不同方法的计算工作量、计算结果和优缺点。 二.程序流程图及图形 1.拉格朗日插值法 开始 x,y,x0 Length (x)==l Ength (y)? n=length (x) i=1:n,l=1。 j=1:i-1&j=i+1:n l=l.*(x0-x(j)/x(i)-x(j) f=f+l*y(i) 结束 否 是 机翼 下轮廓线 2.牛顿插值法 开始 x,y,xi Length(x)==l ength(y)? n=length(x)Y=zeros （n),Y （：1）=y,f=0 a=1:n-1,b=1:n-a,Y(b,a+1)=（Y （b+1，a)-Y(b,a))/(x （b+a)-x(b)) i=1:n,z=1 结束 j=1:i-1,z=z.*(xi-x(j)) f=f+Y(1,i)*z 否 是 3.分段线性插值法 开始 x ，y ，x0 length （x ）==length(y)? k=1:n-1 x(k)<=x0&x0《=x(k+1) temp=x(k)-x(k+1) f=(x0-x(k+1))/temp*y(k)+(x0-x(k))/(-temp)*y(k+1) 结束 否否 是 是 三.matlab 程序及简要的注释（m 文件） 1.拉格朗日插值法 2.牛顿插值法 function f=newdun(x,y,xi) %x 为已知数据点的x 坐标向量 %y 为已知数据点的y 坐标向量 function f=lang(x,y,x0) %x 为已知数据点的x 坐标向量 %y 为已知数据点的y 坐标向量 1电力系统潮流计算 潮流计算是电力系统分析中的一种最基本的计算，它的任务是对给定的运行条件确定系统的运行状态，如母线上的电压(幅值及相角)、网络中的功率分布及功率损耗等。在电力系统规划设计和现有电力系统运行方式的研究中，都需要利用潮流计算来定量地分析比较供电方案或运行方式的合理性.可靠性和经济性。此外，电力系统潮流计算也是计算系统动态稳定和静态稳定的基础。 2节点导纳矩阵的形成 在图1（a ）的简单电力系统中，若略去变压器的励磁功率和线路电容，负荷用阻抗表示，便可以得到一个有5个节点（包括零电位点）和7条支路的等值网络，如图1（b ）所示。将接于节点1和4的电势源和阻抗的串联组合变换成等值的电流源和导纳的并联组合，变得到图1（c ）的等值网络，其中1101I y E =和 4404I y E =分别称为节点1和4的注入电流源。 (a) 2 4? ?4 y (c) 图1 电力系统及其网络 以零电位点作为计算节点电压的参考点，根据基尔霍夫定律，可以写出4个独立节点的电流平衡方程如下： 10112121 12212022323242423323434244234434044()()()()0()()0()()y U y U U I y U U y U y U U y U U y U U y U U y U U y U U y U I ?+-=? -++-+-=? ? -+-=? ?-+-+=? （2-1） 上述方程组经过整理可以写成 1111221211222233244322333344422433444400Y U Y U I Y U Y U Y U Y U Y U Y U Y U Y U Y U Y U I ? + =? +++=? ? ++=?? ++=? （2-2） 式中， 111012 Y y y =+； 2220232412 Y y y y y =+++；332334Y y y =+；44402434Y y y y =++； 122112Y Y y ==-； 233223 Y Y y ==-； 244224 Y Y y ==-； 344334 Y Y y ==-。 一般的，对于有n 个独立节点的网络，可以列写n 个节点方程 11112211211222221122n n n n n n nn n n Y U Y U Y U I Y U Y U Y U I Y U Y U Y U I ? +++=? +++=? ? ? ?++ +=? （2-3） 也可以用矩阵写成 1111121212222212n n n n nn n n U I Y Y Y Y Y Y U I Y Y Y U I ???? ???????? ??????=?????? ?????? ?????????? （2-4） 或缩写为 YU I = （2-5） 矩阵Y 称为节点导纳矩阵。它的对角线元素ii Y 称为节点i 的自导纳，其值等于接于节点i 的所有支路导纳之和。非对角线元素 ij Y 称为节点i 、j 间的互导纳， 它等于直接接于节点i 、j 间的支路导纳的负值。若节点i 、j 间不存在直接支路，则有 ij Y =。由此可知节点导纳矩阵是一个稀疏的对称矩阵。 %均差牛顿插值 function [ f y f0 ] = newton1( X,Y,x0 ) if nargin<3 error('Requires at least three input arguments.'); end if length(X)==length(Y) n=length(X); else error('length must equal') end syms x s=Y(1); l=1.0; y=zeros(n); y(1:n,1)=Y'; for i=2:n for j=2:i y(i,j)=(y(i,j-1)-y(j-1,j-1))/(X(i)-X(j-1)); if i==j l=l*(x-X(i-1)); s=s+y(i,i)*l; end end end f=simple(s); f0=subs(f,x0); function [ f f0 y] = newton2( X,Y,x0 ) if nargin<3 error('Requires at least three input arguments.'); end if length(X)==length(Y) n=length(X); else error('length must equal') end syms x s=Y(1); l=1.0; y=zeros(n) y(1:n,1)=Y'; for i=2:n for j=2:i y(i,j)=(y(i,j-1)-y(i-1,j-1))/(X(i)-X(i-j+1)); if i==j l=l*(x-X(i-1)); s=s+y(i,i)*l; end end end f=simple(s); f0=subs(f,x0); /*利用牛顿-拉夫逊迭代法(极坐标形式)，计算复杂电力系统潮流，具有收敛性好，收敛速度快等优点。所有参数应归算至标幺值下。 /*可计算最大节点数为100，可计算PQ,PV,平衡节点*/ /*可计算非标准变比和平行支路*/ #include 插值MATLAB程序（可以输出多项式）—数值分析 1.拉格朗日多项式逼近 function [C,L,y]=lagran(X,Y) %拉格朗日多项式逼近 w=length(X); L=zeros(w,w); for k=1:w V=1; for j=1:w if k~=j V=conv(V,poly(X(j)))/(X(k)-X(j)); end end L(k,:)=V; end C=Y*L; y=poly2sym(C,'x'); 2.牛顿插值多项式 function [C,D,y]=newpoly(X,Y) %牛顿插值多项式 n=length(X); D=zeros(n,n); D(:,1)=Y'; for j=2:n for k=j:n D(k,j)=(D(k,j-1)-D(k-1,j-1))/(X(k)-X(k-j+1)); end end C=D(n,n); for k=(n-1):-1:1 C=conv(C,poly(X(k))); m=length(C); C(m)=C(m)+D(k,k); end y=poly2sym(C,'x'); 3.切比雪夫逼近 function [C,X,Y]=cheby(fun,n,a,b) %切比雪夫逼近 if nargin==2 a=-1;b=1; end d=pi/(2*n+2); C=zeros(1,n+1); for k=1:n+1 X(k)=cos((2*k-1)*d); end X=(b-a)*X/2+(a+b)/2; x=X; Y=eval(fun); for k=1:n+1 z=(2*k-1)*d; for j=1:n+1 C(j)=C(j)+Y(k)*cos((j-1)*z); end end C=2*C/(n+1); C(1)=C(1)/2;牛顿插值法原理及应用
牛顿插值法matlab程序解析
基于极坐标的牛顿拉夫逊潮流计算
matlab实现数值分析报告插值及积分
牛顿拉夫逊法潮流计算
PQ分解法潮流计算实验
拉格朗日插值、牛顿插值的matlab代码
牛顿——拉夫逊法进行潮流计算
matlab牛顿插值法例题与程序
牛顿插值MATLAB算法
牛顿插值法的MATLAB综合程序
matlab_牛顿插值法_三次样条插值法
matlab计算拉格朗日牛顿及分段线性插值的程序
直角坐标系下牛顿法潮流计算
均差牛顿插值MATLAB,M文件
牛顿-拉夫逊迭代法极坐标潮流计算C语言程序
插值MATLAB程序-数值分析
- MATLAB牛顿插值法例题与程序
- 牛顿插值法matlab代码
- 均差牛顿插值MATLAB,M文件
- matlab 牛顿插值法 三次样条插值法
- matlab-牛顿插值法-三次样条插值法
- MATLAB实现拉格朗日插值
- 牛顿插值MATLAB算法
- matlab(迭代法-牛顿插值)Word版
- matlab迭代法牛顿插值(终审稿)
- Newton插值的matlab实现
- 牛顿插值matlab程序
- matlab(迭代法_牛顿插值)
- MATLAB 牛顿插值法例题与程序
- MATLAB牛顿插值
- 功能强大的牛顿插值matlab程序代码
- matlab牛顿插值法例题与程序
- 插值MATLAB程序-数值分析
- 牛顿插值法的MATLAB综合程序
- matlab 牛顿插值多项式 代码
- matlab实现数值分析报告插值及积分