Question

微分方程数值解

用二阶龙格库塔方法求微分方程Dy/dx=1-(2*x*y)/(1+x*x)', 'y(0)=0',的数值解
我得到的程序是：
MATLAB主程序
function [k,X,Y,fxy,wch,wucha,P]=RK2(funfcn,fun,x0,b,C,y0,h)
x=x0; y=y0;p=128;
n=fix((b-x0)/h);
fxy=zeros(p,1);
wucha=zeros(p,1);
wch=zeros(p,1);
X=zeros(p,1);
Y=zeros(p,length(y));
k=1; X(k)=x; Y(k,:)=y';
% 绘图.
clc,x,h,y
%计算
%fxy=fxy(:);
for k=2:n+1
x=x+h;a2=C(3);b21=C(4);
c1=C(1);
c2=C(2);
x1=x+a2*h;
k1=feval(funfcn,x,y);
y1=y+b21*h*k1;
k2=feval(funfcn,x1,y1);
fxy(k)=feval(fun,x);
y=y+h*(c1*k1+c2*k2);
X(k)=x;
Y(k,:)=y;k=k+1;
plot(X,Y,'mh',X,fxy,'bo')
grid,xlabel('自变量 X'),ylabel('因变量 Y')
legend('用二阶龙格-库塔方法计算dy/dx=f(x,y)，y(x0)=y0在[x0,b]上的数值解','y/dx=f(x,y)，y(x0)=y0的精确解y=f(x)')
end
%计算误差.
for k=2:n+1
wucha(k)=norm(Y(k-1)-Y(k));
wch(k)=norm(fxy(k)-Y(k));
end
X=X(1:k);
Y=Y(1:k,:);
fxy=fxy(1:k,:);
n=1:k;
wucha=wucha(1:k,:);
wch=wch(1:k,:);
P=[n',X,Y,fxy,wch,wucha];
在MATLAB工作窗口输入下面的程序
>>x0=0;b=2;C=[1/4,3/4,2/3,2/3];y0=0;h=1/4;
[k,X,Y,fxy,wch,wucha,P]=RK2(@f1-(2*x*y)/(1+x*x),@(x+1/3*x^3)/(1+x^2),x0,b,C,y0,h)
但得到的是
??? [k,X,Y,fxy,wch,wucha,P]=RK2(@1-(2*x*y)/(1+x*x),@(x+1/3*x^3)/(1+x^2),x0,b,C,y0,h)
|
Error: Missing MATLAB operator.
求高手解答一下

Answer 1

顾名思义，就是方程的数值结果。

微分方程的解，分为解析解和数值解，前者反映的是微分方程的解，可以用一个函数表示；后者同常不能表为初等函数，但是很多问题，我们并不需要解析解，而是能求出一个数值结果就满足了。

举例说，我们希望知道，一个质点从竖直平面内的光滑半圆轨道一端，从静止开始下滑，求质点转过45度经历的时间.这个问题导致一个貌似很简单的一个微分方程:

y'=1/sqrt(sin(x)),即导函数为正选函数平方根的倒数，其解析解不能表示为初等函数形式，但是对于这个问题，我们倒是可以得到任意精确的数值解。

追问

您能告诉我，上面的程序有什么问题吗？我就是在求数值解，可是运行不下去，出现Error: Missing MATLAB operator.的提示是上面问题