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

用matlab编写拉格朗日插值算法的程序
10
[ 标签：matlab,插值算法,程序 ]
用matlab编写拉格朗日插值算法的程序，并以下面给出的函数表为数据基础，在整个插值区间上采用拉格朗日插值法计算f（0.6），写出程序源代码，输出计算结果
x   -2.15     -1.00   0.01   1.02     2.03      3.25
y   17.03   7.24     1.05    2.03    17.06    23.05
匿名 回答:1 人气:6 解决时间:2011-05-24 19:58
满意答案
好评率：83%
做了一个测试，希望有所帮助。
代码：
% 用matlab编写拉格朗日插值算法的程序，并以下面给出的函数表为数据基础，
% 在整个插值区间上采用拉格朗日插值法计算f（0.6），写出程序源代码，输出计算结果
% x   -2.15     -1.00   0.01   1.02     2.03      3.25
% y   17.03   7.24     1.05    2.03    17.06    23.05
function main()
clc;
x = [-2.15     -1.00   0.01   1.02     2.03      3.25];
y = [17.03   7.24     1.05    2.03    17.06    23.05 ];
x0 = 0.6;
f = Language(x,y,x0)
function f = Language(x,y,x0)
%求已知数据点的拉格朗日插值多项式
%已知数据点的x坐标向量: x
%已知数据点的y坐标向量: y
%插值点的x坐标: x0
%求得的拉格朗日插值多项式或在x0处的插值: f
syms t l;
if(length(x) == length(y))
    n = length(x);
else
    disp('x和y的维数不相等！');
    return;   %检错
end
h=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;
end
simplify(h);
if(nargin == 3)
    f = subs (h,'t',x0); %计算插值点的函数值
else
    f=collect(h);
    f = vpa(f,6); %将插值多项式的系数化成6位精度的小数
end
结果：

f =
    0.0201
>>

如何用MATLAB编写的拉格朗日插值算法的程序、二阶龙格-库塔方法的程序和SOR迭代法的程序，要能运行的
∮初夏戀雨￠ 回答:2 人气:29 解决时间:2009-12-08 19:04
满意答案
好评率：100%
拉格朗日
function y=lagrange(x0,y0,x)
n=length(x0);m=length(x);
for i=1:m
z=x(i);
s=0.0;
for k=1:n
p=1.0;
for j=1:n
if j~=k
p=p*(z-x0(j))/(x0(k)-x0(j));
end
end
s=p*y0(k)+s;
end
y(i)=s;
end
SOR迭代法的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;
% 迭代过程
while k<=N
x(1)=(b(1)-A(1,2:n)*x0(2:n)')/A(1,1);
for i=2:n
x(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);
end
if max(abs(x-x0))<=tol
fid = 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;
end
k=k+1;
x0=x;
end
if k==N+1
fid = fopen('SOR_iter_result.txt', 'wt');
fprintf(fid,'\n********用SOR迭代求解线性方程组的输出结果********\n\n');
fprintf(fid,'迭代次数: %d次\n\n',k);
fprintf(fid,'超过最大迭代次数，求解失败！');
fclose(fid);
end
Matlab中龙格-库塔(Runge-Kutta)方法原理及实现龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法。由于此算法精度高，采取措施对误差进行抑制，所以其实现原理也较复杂。该算法是构建在数学支持的基础之上的。龙格库塔方法的理论基础来源于泰勒公式和使用斜率近似表达微分，它在积分区间多预计算出几个点的斜率，然后进行加权平均，用做下一点的依据，从而构造出了精度更高的数值积分计算方法。如果预先求两个点的斜率就是二阶龙格库塔法，如果预先取四个点就是四阶龙格库塔法。一阶常微分方程可以写作：y'=f(x,y),使用差分概念。
(Yn+1-Yn)/h= f(Xn,Yn)推出（近似等于，极限为Yn'）
Yn+1=Yn+h*f(Xn,Yn)
另外根据微分中值定理，存在0<t<1,使得
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;%赋初值，可以是向量，但是要注意维数
for ii=1:n

x(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
调用的子函数以及其调用语句：
function dy=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,[0 15],[1 1 3]);
subplot(121)
plot(T,F)%Matlab自带的ode45函数效果
title('ode45函数效果')
[T1,F1]=runge_kutta1(@test_fun,[1 1 3],0.25,0,15);%测试时改变test_fun的函数维数，别忘记改变初始值的维数
subplot(122)
plot(T1,F1)%自编的龙格库塔函数效果
title('自编的 龙格库塔函数')

用Runge-Kutta方法求一阶微分方程组初值问题的matlab通用程序
2007-11-11 11:43
题目：
用经典Runge-Kutta方法求下列一阶微分方程组的近似解：
y1' = 3y1 + 2y2 − (2x2 + 1)e2x, y1(0) = 1             e2x表示exp（2*x）
y2 '= 4y1 + y2 + (x2 + 2x − 4)e2x, y2(0) = 1
y3 '= 2y1 − y2 − xe3x, y3(0) = 1
y4 '= y1 + x2ex, y4(0) = 1
y5 '= y2 − e2x, y5(0) = 1
其中初值条件y0为一个五维数组，包含了这五个方程在区间[0,1]左端点0的值，并假设N＝10为区间等分数
程序：
function rk(A,x,h,y0)
%A为字符串函数的元胞数组
%x为x轴上的各点
%h为点距
%y0为初值
i=1;
y(i,:)=y0;
m=length(A);
b=x(length(x));
while x(i)<b
a=[x(i),y(i,:)];
for l=1:m
k1(l)=eval(A{l},a);
end
a=[x(i)+h/2,y(i,:)+h/2*k1];
for l=1:l
k2(l)=eval(A{l},a);
end
a=[x(i)+h/2,y(i,:)+h/2*k2];
for l=1:m
k3(l)=eval(A{l},a);
end
a=[x(i)+h,y(i,:)+h*k3];
for l=1:m
k4(l)=eval(A{l},a);
end
y(i+1,:)=y(i,:)+h/6*(k1+2*k2+2*k3+k4);
i=i+1;
end
y
输入及结果：
>> 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 =
[1x38 char] [1x41 char] [1x28 char] [1x21 char] 'a(3)-exp(2*a(1))'
>> rk(A,0:0.1:1,0.1,ones(1,5))
y =
1.0000 1.0000 1.0000 1.0000 1.0000
1.4692 1.1649 1.1312 1.1226 0.9963
2.1246 1.5112 1.3309 1.3031 0.9930
3.0680 2.1507 1.6105 1.5679 1.0078
4.4629 3.2638 1.9873 1.9573 1.0718
6.5725 5.1403 2.4884 2.5335 1.2376
9.8237 8.2476 3.1587 3.3940 1.5930
14.9117 13.3402 4.0745 4.6925 2.2837
22.9721 21.6384 5.3681 6.6746 3.5497
35.8640 35.1212 7.2699 9.7346 5.7843
56.6365 57.0045 10.1818 14.5095 9.6315

matlab 解微分方程组 dx/dt=x+y dy/dt=x-y
2010-10-6 15:56
提问者：刘の鱼 | 浏览次数：409次

推荐答案
2010-10-6 20: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: [1x1 sym]
x: [1x1 sym]

>> disp(ans.x)
C1/exp(2^(1/2)*t) + C2*exp(2^(1/2)*t) - (2^(1/2)*C1)/exp(2^(1/2)*t) + 2^(1/2)*C2*exp(2^(1/2)*t)

>> disp(ans.y)
C1/exp(2^(1/2)*t) + C2*exp(2^(1/2)*t)
参考资料：