y1(1)=2;y2(1)=3;x(1)=0; for i=1:100;
K1=f(x(i),y1(i),y2(i)); L1=g(x(i),y1(i),y2(i));
K2=f(x(i)+0.5*h,y1(i)+0.5*h*K1,y2(i)+0.5*h*L1); L2=g(x(i)+0.5*h,y1(i)+0.5*h*K1,y2(i)+0.5*h*L1); K3=f(x(i)+0.5*h,y1(i)+0.5*h*K2,y2(i)+0.5*h*L2); L3=g(x(i)+0.5*h,y1(i)+0.5*h*K2,y2(i)+0.5*h*L2); K4=f(x(i)+h,y1(i)+h*K3,y2(i)+h*L3); L4=g(x(i)+h,y1(i)+h*K3,y2(i)+h*L3); x(i+1)=x(i)+h;
y1(i+1)=y1(i)+h*(1/6)*(K1+2*K2+2*K3+K4); y2(i+1)=y2(i)+h*(1/6)*(L1+2*L2+2*L3+L4); end
disp('[x y1 y2]') [y1 ;y2]'
plot(x,y1,'r-',x,y2,'b--','LineWidth',4) legend('y1','y2');
title('R--K四阶龙格库塔算法下方程组的解'); ylabel('y1曲线 y2曲线')
xlabel('x')
7.用有限差分法求解边值问题(h=0.1):
y (1 x2)y 0
.
y( 1) y(1) 1
clc;clear;
% %有限差分法求微分方程
%*********************************************** % -Y''+q(x)y=f(x)
% y(a)=m, 左端点函数值 % y(b)=n,右端点函数值
%以上是说明部分 常微分方程的形式
%************************************************
%%此题求解的是方程 y''-(1+x^2)*y=0 y(-1)=y(1)=1;