[转帖]同伦算法的Matlab程序

suffer

%% N 为积分部数，h 为积分步长，其值视计算精度而定

%% x0 为初始值，可任意给定

%% Equfun(x)为非线性方程组表达式（列向量），

%% Jacobfun(x)为雅可比矩阵

N=100;
h=1/N;
x0=[0 0 0]';
x=x0;

%format long
f=Equfun(x);
b=-h*f;

for i=1:N
A=Jacobfun(x);
k1=inv(A)*b;
A=Jacobfun(x+0.5*k1);
k2=inv(A)*b;
A=Jacobfun(x+0.5*k2);
k3=inv(A)*b;
A=Jacobfun(x+0.5*k3);
k4=inv(A)*b;
x=x+(k1+2*k2+2*k3+k4)/6;
end

disp('The Solution is:')
disp('x=');disp(x);

%%%%%%%%%%%%%%%%%%%%%%%%%

%%%下面是非线性方程组及其雅可比行列式

%%%%%%%%%%%%%%%%%%%%%%%%%

function y=Equfun(x)

%nonliner functions

y=[3*x(1)-cos(x(2)*x(3))-1/2;
x(1)^2-81*(x(2)+0.1)^2+sin(x(3))+1.06;
exp(-x(1)*x(2))+20*x(3)+(10*pi-3)/3];

%%%%%%%%%%%%%%%%%%%%%%%%%

function y=Jacobfun(x)

%Jacobi function
%J=zeros(length(x));

y(1,1)=3;
y(1,2)=x(3)*sin(x(2)*x(3));
y(1,3)=x(2)*sin(x(2)*x(3));
y(2,1)=2*x(1);
y(2,2)=-162*(x(2)+0.1);
y(2,3)=cos(x(3));
y(3,1)=-x(2)*exp(-x(1)*x(2));
y(3,2)=-x(1)*exp(-x(1)*x(2));
y(3,3)=20;

siyanger

谢谢分享！同伦算法也可解常微分方程组吧？

suffer

原帖由 siyanger 于 2006-6-30 07:22 发表
谢谢分享！同伦算法也可解常微分方程组吧？

好像不能直接用于常微分方程求解

21172485

程序验证过了,结果不正确

yuyugongwu

哪有这么简单啊，给大家推荐一本书：非线性数值分析，上面有详细的算法，很好写。

hopeasy

楼上的说的太泛了，《非线性数值分析》的书名字很多，要具体到作者，出版社或者出版年份