马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
求助方程组的边值问题
function [T,Z]=rks41(F,a,b,Za,M)
%Input -F is the system input as a string 'F'
% -a and b are the end points of the interval
% -Za=[x(a) y(a)] are the initial conditions
% -M is the number of steps
%Output -T is the vector of steps
% -Z=[x1(t) ...xn(t)]; where xk(t) is the approximation
% to the kth dependent variable
h=(b-a)/M;
T=zeros(1,M+1);
Z=zeros(M+1,length(Za));
T=a:h:b;
Z(1,:)=Za;
for j=1:M
k1=h*feval(F,T(j),Z(j,:));
k2=h*feval(F,T(j)+h/2,Z(j,:)+k1/2);
k3=h*feval(F,T(j)+h/2,Z(j,:)+k2/2);
k4=h*feval(F,T(j)+h,Z(j,:)+k3);
Z(j+1,:)=Z(j,:)+(k1+2*k2+2*k3+k4)/6;
end
function Z=F1(x,Z)
y1=Z(1);
y2=Z(2);
y3=Z(3);
y4=Z(4);
Z=[y2,-2*y1*y3-2*3*cos(4*x)*y1+2*1.2*y1,y4,0.5*(y3-y1^2)];
%Z=[y2;-2*y1+3*y2;y4;-2*y3+3*y4];
function Z=F2(x,Z)
y1=Z(1);
y2=Z(2);
y3=Z(3);
y4=Z(4);
Z=[y2,-2*y1*y3-2*3*cos(4*x)*y1+2*1.2*y1,y4,0.5*(y3)];
%Z=[y2;x-2*y1+3*y2;y4;x-2*y3+3*y4];
function L=linsht1(F1,F2,a,b,alpha1,alpha2,beta1,beta2,M)
%线性打靶法 求解非线性方程组的边值问题
%Input -F1 and F2 are the systems of first-order equations
% representing the I.V.P.'s (1) and (2), respectively;
% 设w(t),n(t)为边值问题
% w''(t)=-2*w(t)*n(t)-2*3*cos(4*t)*w(t)+2*1.2*w(t), w(a)=0, w(b)=0 (1)
% n''(t)=0.5*(n(t)-w(t)^2), n(a)=0, n(b)=0 (2)
%
% input as strings 'F1','F2'
% - a and b are the end points of the interval
% -alpha=x(a) and beta=x(b); boundary conditions
% -M is the number of steps
% Output -L=[T' X]; where T' is the (M+1)x1 vector of ordinates
%Solve the system F1
Za=[alpha1,0,alpha2,0];
[T,Z]=rks41(F1,a,b,Za,M);
U1=Z(:,1);
U2=Z(:,3);
%Solve the sysstem F2
Za=[0,1,0,1];
[T,Z]=rks41(F2,a,b,Za,M);
V1=Z(:,1);
V2=Z(:,3);
%Calculate the solution to the boundary value problem
X1=U1+(beta1-U1(M+1))*V1/V1(M+1);
X2=U2+(beta2-U2(M+1))*V2/V2(M+1);
%L1=[T' U1 U2];
L=[T' X1 X2];
程序运行无错误,可是结果不正确,大家帮忙看看有什么问题,如有做过类似程序,能否提供参考(szwstar@163.com),谢谢 |