用matlab编写拉格朗日插值算法的程序

用matlab编写拉格朗日插值算法的程序10[标签：matlab,插值算法,程序]用matlab编写拉格朗日插值算法的程序，并以下面给出的函数表为数据基础，在整个插值区间上采用拉格朗日插值法计算f（0.6），写出程序源代码，输出计算结果x-2.15-1.000.011.022.033.25y17.037.241.052.0317.0623.05匿名回答:1人气:6解决时间:2011-05-2419:58满意答案好评率：83%做了一个测试，希望有所帮助。代码：%用matlab编写拉格朗日插值算法的程序，并以下面给出的函数表为数据基础，%在整个插值区间上采用拉格朗日插值法计算f（0.6），写出程序源代码，输出计算结果%x-2.15-1.000.011.022.033.25%y17.037.241.052.0317.0623.05functionmain()clc;x=[-2.15-1.000.011.022.033.25];y=[17.037.241.052.0317.0623.05];x0=0.6;f=Language(x,y,x0)functionf=Language(x,y,x0)%求已知数据点的拉格朗日插值多项式%已知数据点的x坐标向量:x%已知数据点的y坐标向量:y%插值点的x坐标:x0%求得的拉格朗日插值多项式或在x0处的插值:fsymstl;if(length(x)==length(y))n=length(x);elsedisp('x和y的维数不相等！');return;%检错endh=sym(0);for(i=1:n)l=sym(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;h=h+l;endsimplify(h);if(nargin==3)f=subs(h,'t',x0);%计算插值点的函数值elsef=collect(h);f=vpa(f,6);%将插值多项式的系数化成6位精度的小数end结果：f=0.0201如何用MATLAB编写的拉格朗日插值算法的程序、二阶龙格-库塔方法的程序和SOR迭代法的程序，要能运行的∮初夏戀雨￠回答:2人气:29解决时间:2009-12-0819:04满意答案好评率：100%拉格朗日functiony=lagrange(x0,y0,x)n=length(x0);m=length(x);fori=1:mz=x(i);s=0.0;fork=1:np=1.0;forj=1:nifj~=kp=p*(z-x0(j))/(x0(k)-x0(j));endends=p*y0(k)+s;endy(i)=s;endSOR迭代法的Matlab程序function[x]=SOR_iterative(A,b)%用SOR迭代求解线性方程组,矩阵A是方阵x0=zeros(1,length(b));%赋初值tol=10^(-2);%给定误差界N=1000;%给定最大迭代次数[n,n]=size(A);%确定矩阵A的阶w=1;%给定松弛因子k=1;%迭代过程whilek=Nx(1)=(b(1)-A(1,2:n)*x0(2:n)')/A(1,1);fori=2:nx(i)=(1-w)*x0(i)+w*(b(i)-A(i,1:i-1)*x(1:i-1)'-A(i,i+1:n)*x0(i+1:n)')/A(i,i);endifmax(abs(x-x0))=tolfid=fopen('SOR_iter_result.txt','wt');fprintf(fid,'\n********用SOR迭代求解线性方程组的输出结果********\n\n');fprintf(fid,'迭代次数:%d次\n\n',k);fprintf(fid,'x的值\n\n');fprintf(fid,'%12.8f\n',x);break;endk=k+1;x0=x;endifk==N+1fid=fopen('SOR_iter_result.txt','wt');fprintf(fid,'\n********用SOR迭代求解线性方程组的输出结果********\n\n');fprintf(fid,'迭代次数:%d次\n\n',k);fprintf(fid,'超过最大迭代次数，求解失败！');fclose(fid);endMatlab中龙格-库塔(Runge-Kutta)方法原理及实现龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法。由于此算法精度高，采取措施对误差进行抑制，所以其实现原理也较复杂。该算法是构建在数学支持的基础之上的。龙格库塔方法的理论基础来源于泰勒公式和使用斜率近似表达微分，它在积分区间多预计算出几个点的斜率，然后进行加权平均，用做下一点的依据，从而构造出了精度更高的数值积分计算方法。如果预先求两个点的斜率就是二阶龙格库塔法，如果预先取四个点就是四阶龙格库塔法。一阶常微分方程可以写作：y'=f(x,y),使用差分概念。(Yn+1-Yn)/h=f(Xn,Yn)推出（近似等于，极限为Yn'）Yn+1=Yn+h*f(Xn,Yn)另外根据微分中值定理，存在0t1,使得Yn+1=Yn+h*f(Xn+th,Y(Xn+th))这里K＝f(Xn+th,Y(Xn+th))称为平均斜率，龙格库塔方法就是求得K的一种算法。利用这样的原理，经过复杂的数学推导（过于繁琐省略），可以得出截断误差为O(h^5)的四阶龙格库塔公式：K1＝f(Xn,Yn);K2=f(Xn+h/2,Yn+(h/2)*K1);K3=f(Xn+h/2,Yn+(h/2)*K2);K4=f(Xn+h,Yn+h*K3);Yn+1=Yn+h*(K1+2K2+2K3+K4)*(1/6);所以，为了更好更准确地把握时间关系，应自己在理解龙格库塔原理的基础上，编写定步长的龙格库塔函数，经过学习其原理，已经完成了一维的龙格库塔函数。仔细思考之后，发现其实如果是需要解多个微分方程组，可以想象成多个微分方程并行进行求解，时间，步长都是共同的，首先把预定的初始值给每个微分方程的第一步，然后每走一步，对多个微分方程共同求解。想通之后发现，整个过程其实很直观，只是不停的逼近计算罢了。编写的定步长的龙格库塔计算函数：function[x,y]=runge_kutta1(ufunc,y0,h,a,b)%参数表顺序依次是微分方程组的函数名称，初始值向量，步长，时间起点，时间终点（参数形式参考了ode45函数）n=floor((b-a)/h);%求步数x(1)=a;%时间起点y(:,1)=y0;%赋初值，可以是向量，但是要注意维数forii=1:nx(ii+1)=x(ii)+h;k1=ufunc(x(ii),y(:,ii));k2=ufunc(x(ii)+h/2,y(:,ii)+h*k1/2);k3=ufunc(x(ii)+h/2,y(:,ii)+h*k2/2);k4=ufunc(x(ii)+h,y(:,ii)+h*k3);y(:,ii+1)=y(:,ii)+h*(k1+2*k2+2*k3+k4)/6;%按照龙格库塔方法进行数值求解end调用的子函数以及其调用语句：functiondy=test_fun(x,y)dy=zeros(3,1);%初始化列向量dy(1)=y(2)*y(3);dy(2)=-y(1)+y(3);dy(3)=-0.51*y(1)*y(2);对该微分方程组用ode45和自编的龙格库塔函数进行比较，调用如下：[T,F]=ode45(@test_fun,[015],[113]);subplot(121)plot(T,F)%Matlab自带的ode45函数效果title('ode45函数效果')[T1,F1]=runge_kutta1(@test_fun,[113],0.25,0,15);%测试时改变test_fun的函数维数，别忘记改变初始值的维数subplot(122)plot(T1,F1)%自编的龙格库塔函数效果title('自编的龙格库塔函数')用Runge-Kutta方法求一阶微分方程组初值问题的matlab通用程序2007-11-1111:43题目：用经典Runge-Kutta方法求下列一阶微分方程组的近似解：y1'=3y1+2y2−(2x2+1)e2x,y1(0)=1e2x表示exp（2*x）y2'=4y1+y2+(x2+2x−4)e2x,y2(0)=1y3'=2y1−y2−xe3x,y3(0)=1y4'=y1+x2ex,y4(0)=1y5'=y2−e2x,y5(0)=1其中初值条件y0为一个五维数组，包含了这五个方程在区间[0,1]左端点0的值，并假设N＝10为区间等分数程序：functionrk(A,x,h,y0)%A为字符串函数的元胞数组%x为x轴上的各点%h为点距%y0为初值i=1;y(i,:)=y0;m=length(A);b=x(length(x));whilex(i)ba=[x(i),y(i,:)];forl=1:mk1(l)=eval(A{l},a);enda=[x(i)+h/2,y(i,:)+h/2*k1];forl=1:lk2(l)=eval(A{l},a);enda=[x(i)+h/2,y(i,:)+h/2*k2];forl=1:mk3(l)=eval(A{l},a);enda=[x(i)+h,y(i,:)+h*k3];forl=1:mk4(l)=eval(A{l},a);endy(i+1,:)=y(i,:)+h/6*(k1+2*k2+2*k3+k4);i=i+1;endy输入及结果：y1='3*a(2)+2*a(3)-(2*a(1)^2+1)*exp(2*a(1))';解释：a(1)为x，a(2)、a(3)……为y1，y2……y2='4*a(2)+a(3)+(a(1)^2+2*a(1)-4)*exp(2*a(1))';y3='2*a(2)-a(3)-a(1)*exp(3*a(1))';y4='a(2)+a(1)^2*exp(a(1))';y5='a(3)-exp(2*a(1))';A=[{y1},{y2},{y3},{y4},{y5}]A=[1x38char][1x41char][1x28char][1x21char]'a(3)-exp(2*a(1))'rk(A,0:0.1:1,0.1,ones(1,5))y=1.00001.00001.00001.00001.00001.46921.16491.13121.12260.99632.12461.51121.33091.30310.99303.06802.15071.61051.56791.00784.46293.26381.98731.95731.07186.57255.14032.48842.53351.23769.82378.24763.15873.39401.593014.911713.34024.07454.69252.283722.972121.63845.36816.67463.549735.864035.12127.26999.73465.784356.636557.004510.181814.50959.6315matlab解微分方程组dx/dt=x+ydy/dt=x-y2010-10-615:56提问者：刘の鱼|浏览次数：409次推荐答案2010-10-620:50不知道解得对不对程序：dsolve('Dx=x+y','Dy=x-y','t')解得：x=C1*exp(2^(1/2)*t)+C2*exp(-2^(1/2)*t)y=C1*2^(1/2)*exp(2^(1/2)*t)-C2*2^(1/2)*exp(-2^(1/2)*t)-C1*exp(2^(1/2)*t)-C2*exp(-2^(1/2)*t)dsolve('Dx=x+y','Dy=x-y')ans=y:[1x1sym]x