%书籍：常用数值算法及其matlab实现%第10章 常微分方程初值问题的数值解法,例10.14使用%四阶龙格库塔方法function [t,z] = rk4symeq(fun, t0, tf, Za, h)%fun:微分方程的右表达式%t0, tn为区间%Za为初值,是列向量M = floor(tf-t0)/h ;%离散点的个数M+1if t0 >= tfprintf('左端点必须小于右端点');return;endN = length(Za); %获得变量个数,Nz = zeros(M+1, N);t = zeros(M +1, 1);t =[t0 : h :tf]';z(1,:) = Za'; %假设Za为列向量，与微分方程中的变量方向统一，变成行向量for i = 1:MK1 = feval(fun, t(i) , z(i,:));%K是行向量K2 = feval(fun, t(i)+1/2*h ,z(i,:)+1/2* h*K1);K3 = feval(fun, t(i)+1/2*h ,z(i,:)+1/2* h*K2);K4 = feval(fun, t(i)+ h ,z(i,:)+ h*K3); z(i+1,:) = z(i,:) +h/6 *(K1 + 2*K2 + 2*K3 + K4);end以下为求解的方程组%书籍：常用数值算法及其matlab实现%第10章 常微分方程初值问题的数值解法%四阶龙格库塔，例10.16function s = exa10_16(t,z)%z是个向量,1*3%输出s也是向量,1*3s = zeros(1,2);dy1 = z(2) ;dy2 = -z(1) +(1-z(1)^2)*z(2);s = [dy1 dy2 ];%输出s为行向量主函数clear all;clc;close all;%书籍：常用数值算法及其matlab实现%第10章 常微分方程组初值问题的数值解法%四阶龙格库塔，例10.16%函数原型 function [t,z] = rk4symeq(fun, t0, tn, Za, h)format longt0 = 0; tf = 15;Za = [-3 ; -0.1];%x初值h1 = 0.01; [t1,z1] = rk4symeq(@exa10_16, t0, tf , Za, h1);figure(1)plot(t1,z1(:,1),'b',t1,z1(:,2), 'r*')legend('y1','y2')figure(2)plot(z1(:,1),z1(:,2));xlabel('x');ylabel('y');
运行结果：
