




版权说明:本文档由用户提供并上传,收益归属内容提供方,若内容存在侵权,请进行举报或认领
文档简介
1、MATLAB程序设计实践课程考核一、编程实现“四阶龙格库塔(R-K)方法求常微分方程”,并举一例应用之。【实例】采用龙格-库塔法求微分方程:1、算法说明:在龙格-库塔法中,四阶龙格-库塔法的局部截断误差约为o(h5),被广泛应用于解微分方程的初值问题。其算法公式为:其中:2、流程图:2.1、四阶龙格库塔(R-K)方法流程图:输入待求微分方程、求解的自变量范围、初值以及求解范围内的取点数等。确定求解范围内的步长k = 取点数?否求解:求解并输出:是结束算法2.2、实例求解流程图:输入求解的自变量范围求出待求简单微分方程的真值解用MATLAB自带函数ode23求解待求微分方程结束用自编函数四阶龙格
2、库塔(R-K)方法求解待求微分方程开始3、源程序代码3.1、四阶龙格库塔(R-K)方法源程序:function x,y = MyRunge_Kutta(fun,x0,xt,y0,PointNum,varargin)%Runge-Kutta 方法解微分方程形为 y'(t)=f(x,y(x)%此程序可解高阶的微分方程。只要将其形式写为上述微分方程的向量形式%函数 f(x,y): fun%自变量的初值和终值:x0, xt%y0表示函数在x0处的值,输入初值为列向量形式%自变量在x0,xt上取的点数:PointNum%varargin为可输入项,可传适当参数给函数f(x,y)%x:所取的点的x
3、值%y:对应点上的函数值 if nargin<4 | PointNum<=0 PointNum=100;endif nargin<3 y0=0;endy(1,:)=y0(:)' %初值存为行向量形式h=(xt-x0)/(PointNum-1); %计算步长 x=x0+0:(PointNum-1)'*h; %得x向量值for k=1:(PointNum) %迭代计算 f1=h*feval(fun,x(k),y(k,:),varargin:); f1=f1(:)' %得公式k1 f2=h*feval(fun,x(k)+h/2,y(k,:)+f1/2,var
4、argin:); f2=f2(:)' %得公式k2 f3=h*feval(fun,x(k)+h/2,y(k,:)+f2/2,varargin:); f3=f3(:)' %得公式k3 f4=h*feval(fun,x(k)+h,y(k,:)+f3,varargin:); f4=f4(:)' %得公式k4 y(k+1,:)=y(k,:)+(f1+2*(f2+f3)+f4)/6; %得y(n+1)end3.2、实例求解源程序:%运行四阶R-K法clear, clc %清除内存中的变量x0=0;xt=2;Num=100;h=(xt-x0)/(Num-1);x=x0+0:Num*
5、h;a=1;yt=1-exp(-a*x); %真值解fun=inline('-y+1','x','y'); %用inline构造函数f(x,y)y0=0; %设定函数初值PointNum=5; %设定取点数x1,y1=ode23(fun,0,2,0);xr,yr=MyRunge_Kutta(fun,x0,xt,y0,PointNum);MyRunge_Kutta_x=xr'MyRunge_Kutta_y=yr'plot(x,yt,'k',x1,y1,'b-',xr,yr,'r-')l
6、egend('真值','ode23','Rung-Kutta法解',2)hold onplot(x1,y1,'bo',xr,yr,'r*')4、程序运行结果:MyRunge_Kutta_x = 0 0.5000 1.0000 1.5000 2.0000MyRunge_Kutta_y = 0 0.3932 0.6318 0.7766 0.8645二、编程解决以下科学计算问题:(一)例7-2-4 材料力学复杂应力状态的分析Moore圆。1、程序说明:利用平面应力状态下斜截面应力的一般公式,画出任意平面应力状态下的应力圆
7、(Moore圆),求出相应平面应力状态下的主应力(、),并求出该应力状态下任意方位角的斜截面上的应力、。2、程序流程图:开始输入待求应力状态的参数画出应力圆求某一方向角截面上的应力?输入方向角求出相应正应力、切应力是否得出该应力状态下的主应力求出主应力平面方向角结束3、程序代码:clear;clc;Sx=input('Sigma_x(MPa)='); %输入该应力状态下的各应力值Sy=input('Sigma_y(MPa)=');Txy=input('Tau_xy(MPa)=');a=linspace(0,pi,100); %等分半圆周角Sa=(
8、Sx+Sy)/2;Sd=(Sx-Sy)/2;Sigma=Sa+Sd*cos(2*a)-Txy*sin(2*a); %应力圆一般方程Tau=Sd*sin(2*a)+Txy*cos(2*a);plot(Sigma,Tau,Sx,Txy,'.r',Sy,-Txy,'.r'); %画出应力圆,标出该应力状态下各应力参数line(Sx,Sy,Txy,-Txy);axis equal; %控制各坐标轴的分度使其相等使应力圆变圆title('应力圆');xlabel('正应力(MPa)');ylabel('剪应力(MPa)');
9、text(Sx,Txy,'A');text(Sy,-Txy,'B');Smax=max(Sigma);Smin=min(Sigma);Tmax=max(Tau);Tmin=min(Tau);b=axis; %提取坐标轴边界ps_array.Color='k' %控制坐标轴颜色为黑色line(b(1),b(2),0,0,ps_array); %调整坐标轴line(0,0,b(3),b(4),ps_array); b=axis; %提取坐标轴边界line(b(1),b(2),0,0,ps_array); %画出x坐标轴hold onplot(Sa,0
10、,'.r') %标出圆心text(Sa,0,'O');plot(Smax,0,'.r',Smin,0,'.r',Sa,Tmax,'.r',Sa,Tmin,'.r') %标出最大、最小拉应力、切应力点text(Smax,0,'C');text(Smin,0,'D');text(Sa,Tmax,'E');text(Sa,Tmin,'F');%-此部分求某一斜截面上的应力-t=input('若需求某一截面上的应力,请输入1;若不求应力
11、,请输入0:');while t=0 alpha=input('给出斜截面方向角:alpha=(角度):'); sigma=Sa+Sd*cos(2*(alpha/180*pi)-Txy*sin(2*(alpha/180*pi) tau=Sd*sin(2*(alpha/180*pi)+Txy*cos(2*(alpha/180*pi) plot(sigma,tau,'or') t=input('若还需求其他截面上的应力,请输入1;若要退出,请输入0:');endhold off%-此部分求该应力状态下的主应力-Sigma_Max=SmaxSi
12、gma_Min=SminTau_Max=TmaxTau_Min=TminSigma1=Smax %得出主应力Sigma3=Sminl=Sx-Sa;h=Txy;ratio=abs(h/l); %求主应力平面方向角'主应力平面方向角(角度):'alpha_0=atan(ratio)/2*180/pi4、程序运行结果:(以为例)Sigma_x(MPa)=100Sigma_y(MPa)=30Tau_xy(MPa)=-20若需求某一截面上的应力,请输入1;若不求应力,请输入0:1给出斜截面方向角:alpha=(角度):30sigma = 99.8205tau = 20.3109若还需求其
13、他截面上的应力,请输入1;若要退出,请输入0:0Sigma_Max = 105.3087Sigma_Min = 24.6970Tau_Max = 40.3109Tau_Min = -40.2963Sigma1 = 105.3087Sigma3 = 24.6970ans =主应力平面方向角(角度):alpha_0 = 14.8724(二)实验5(椭圆的交点) 两个椭圆可能具有0 4个交点,求下列两个椭圆的所有交点坐标:(1) ; (2) 1、算法说明:此题相当于求两个二元二次方程组的解,故为便于清楚地显示出两椭圆的相对位置,用ezplot函数把两个椭圆画在同一个坐标系中,然后利用fsolve函数解方程组得到两椭圆的交点即方程组的解。2、程序流程图:开始依照所给方程画出两个椭圆用符号解法求两椭圆的交点将符号解转换成数值解结束3、程序代码:clear; clc;ezplot('(x-2)2+(y-3+2*x)2-5',-1,5, -8,8); %画第一个椭圆hold onezplot('2*(x-3)2+(y/3)2-4',-1,5, -8,8); %画第二个椭圆grid on; %显示网格hold
温馨提示
- 1. 本站所有资源如无特殊说明,都需要本地电脑安装OFFICE2007和PDF阅读器。图纸软件为CAD,CAXA,PROE,UG,SolidWorks等.压缩文件请下载最新的WinRAR软件解压。
- 2. 本站的文档不包含任何第三方提供的附件图纸等,如果需要附件,请联系上传者。文件的所有权益归上传用户所有。
- 3. 本站RAR压缩包中若带图纸,网页内容里面会有图纸预览,若没有图纸预览就没有图纸。
- 4. 未经权益所有人同意不得将文件中的内容挪作商业或盈利用途。
- 5. 人人文库网仅提供信息存储空间,仅对用户上传内容的表现方式做保护处理,对用户上传分享的文档内容本身不做任何修改或编辑,并不能对任何下载内容负责。
- 6. 下载文件中如有侵权或不适当内容,请与我们联系,我们立即纠正。
- 7. 本站不保证下载资源的准确性、安全性和完整性, 同时也不承担用户因使用这些下载资源对自己和他人造成任何形式的伤害或损失。
最新文档
- 通辽市人民医院医疗行业宏观环境PEST分析理解试题
- 鹤岗市人民医院肝胆外科副主任医师晋升主任医师考核
- 巴彦淖尔市人民医院耳鼻喉头颈外科住院医师规范化培训考核
- 重庆市中医院皮瓣移植护理考核
- 临汾市人民医院水质监测专项考核
- 2025妇幼保健院护理管理体系规划考核
- 2025年造纸厂建设项目可行性分析报告
- 七台河市人民医院介入质量控制考核
- 重庆市人民医院细胞学快速诊断考核
- 中国食品检测试剂项目商业计划书
- GB/T 41782.4-2024物联网系统互操作性第4部分:语法互操作性
- 托班自主活动教案
- 视频新媒体制作服务方案
- 中华民族共同体概论课件专家版2第二讲 树立正确的中华民族历史观
- 宫腔镜检查的个案护理
- 大学生创新创业基础(创新创业课程)全套教学课件
- 初中英语短语大全to-do-结构
- 全套ISO45001职业健康安全管理体系文件(手册及程序文件)
- 敦煌文献研究与敦煌学
- 各种膏药的配方
- 笛卡尔环线性化技术的基本原理
评论
0/150
提交评论