MapleSim in Engineering Education

Maple's dsolve numeric can solve delay ODEs and DAEs as of Maple 18. However, if I am not wrong, it cannot solve delay equations with a time dependent history. In this post I show two examples.

Example 1:

y1(t) and y2(t) with time dependent history. Use of piecewise helps this problem to be solved efficiently. Hopefully Maple will add history soon in its capability.

Example 2: 

This is a very a complicated stiff problem from immunology. As of now, I believe only Maple can solve this (other than RADAR5 from Prof. Hairer). Details and plots are posted in the attached code.

 

Let me know if any one has a delay problem that needs to be solved. I have tested many delay problems in Maple (they work fine). The attached examples required addtional tweaking, hence the post.

 

I want to take this opportunity to congratulate and thank Maple's dsolve numeric/delay solvers for their fantastic job. Maple is world leader not because of example1, but because of its ability to solve example 2.

 

 

restart;

 This code is written by Dayaram Sonawane and Venkat R. Subramnian, University of Washington. You will need Maple 18 or later for this. For those who are wanting to solve these problems in earlier versions, I can help them by offering a procedure based approach (less efficient).

Example1 The first example solved is a state dependent delay problem (http://www.mathworks.com/help/matlab/math/state-dependent-delay-problem.html).

 

eq1:= diff(y1(t),t)=y2(t);

eq1 := diff(y1(t), t) = y2(t)

(1)

eq2:=diff(y2(t),t)=-y2(exp(1-y2(t)))*y2(t)^2*exp(1-y2(t));

eq2 := diff(y2(t), t) = -y2(exp(1-y2(t)))*y2(t)^2*exp(1-y2(t))

(2)

 Both y1(t) and y2(t) have time dependent history (y1(t)=log(t) and y2(t)=1/t, t<-0.1). If I am not mistaken one cannot solve this directly using Maple's dsolve numeric command. However, a simple trick can be used to redefine the equations for y1(t) and y2(t) as below

eq3:=diff(y1(t),t)=piecewise(t<=0.1,1/t,y2(t));

eq3 := diff(y1(t), t) = piecewise(t <= .1, 1/t, y2(t))

(3)

eq4:=diff(y2(t),t)=piecewise(t<=0.1,-1/t^2,-y2(exp(1-y2(t)))*y2(t)^2*exp(1-y2(t)));

eq4 := diff(y2(t), t) = piecewise(t <= .1, -1/t^2, -y2(exp(1-y2(t)))*y2(t)^2*exp(1-y2(t)))

(4)

 The problem is solved from a small number close to t = 0 (1e-4) to make Maple's dsolve numeric remember the history till t = 0.1

epsilon:=1e-4;

epsilon := 0.1e-3

(5)

sol:=dsolve({eq3,eq4,y1(epsilon)=log(epsilon),y2(epsilon)=1/epsilon},type=numeric,delaymax=5):

with(plots):

odeplot(sol,[t,y1(t)],0.1..5,thickness=3,axes=boxed);

 

odeplot(sol,[t,y2(t)],0.1..5,thickness=3,axes=boxed);

 

sol(5.0);log(5.0);1/5.0;

[t = 5.0, y1(t) = 1.60942323180838, y2(t) = .199998786891688]

1.609437912

.2000000000

(6)

Tweaking the tolerances and epsilon will get the solution even more closer to the expected answers.

 

 

 Example 2

 The next problem discussed is very stiff, complicated and as of today, according Professor Hairer (one of the world's leading authority in numerical solutions of ODEs, DAEs), cannot be solved by any other code other than his RADAR (5th order implicit Runge Kutta modified for delay equations, Guglielmi N. and Hairer E. (2001) Implementing Radau IIa methods for stiff delay differential equations. Computing 67:1-12). This problem requires very stringent tolerances. For more information read, http://www.scholarpedia.org/article/Stiff_delay_equations. I can safely say that Maple can boast that it can solve this delay differential equation by using a switch function (instead of Heaviside/picecewise function). Code is attached below and results are compared with the output from RADAR code.  Note that dsolve/numeric is probably taking more time steps compared to RADAR, but the fact that Maple's dsolve numeric solved this model (which cannot be solved in Mathematica or MATLAB[needs confirmation for MATLAB]) should make Maple's code writers proud. It is very likely that we will be trying to submit an educational/research article on this topic/example soon to a journal. For some weird reasons, stiff=true gives slightly inaccurate results.

restart:

 

radar5data:=readdata("C:\\Users\\Venkat16core-office\\Google Drive\\waltmanproblem\\sol.txt",[string,string,float,string,string,float,float,float,float,float,float]):

nops(radar5data);

1059

(7)

radar5data[1059];

["X", "=", 300.000000, "Y", "=", 0.6154486288e-15, 0.3377120916e-6, 0.4221403310e-6, 0.2142554563e-5, 299.9999999, 299.6430338]

(8)

eq[1]:=diff(y[1](t),t)=-r*y[1](t)*y[2](t)-s*y[1](t)*y[4](t);

eq[1] := diff(y[1](t), t) = -r*y[1](t)*y[2](t)-s*y[1](t)*y[4](t)

(9)

eq[2]:=diff(y[2](t),t)=-r*y[1](t)*y[2](t)+alpha*r*y[1](y[5](t))*y[2](y[5](t))*H1;#Heaviside(t-35);

eq[2] := diff(y[2](t), t) = -r*y[1](t)*y[2](t)+alpha*r*y[1](y[5](t))*y[2](y[5](t))*H1

(10)

eq[3]:=diff(y[3](t),t)=r*y[1](t)*y[2](t);

eq[3] := diff(y[3](t), t) = r*y[1](t)*y[2](t)

(11)

eq[4]:=diff(y[4](t),t)=-s*y[1](t)*y[4](t)-gamma1*y[4](t)+beta*r*y[1](y[6](t))*y[2](y[6](t))*H2;#Heaviside(t-197);

eq[4] := diff(y[4](t), t) = -s*y[1](t)*y[4](t)-gamma1*y[4](t)+beta*r*y[1](y[6](t))*y[2](y[6](t))*H2

(12)

eq[5]:=diff(y[5](t),t)=H1*(y[1](t)*y[2](t)+y[3](t))/(y[1](y[5](t))*y[2](y[5](t))+y[3](y[5](t)));#eq[7]:=y[7](t)=HH(t);

eq[5] := diff(y[5](t), t) = H1*(y[1](t)*y[2](t)+y[3](t))/(y[1](y[5](t))*y[2](y[5](t))+y[3](y[5](t)))

(13)

eq[6]:=diff(y[6](t),t)=H2*(10.^(-12)*0+y[2](t)+y[3](t))/(10.^(-12)*0+y[2](y[6](t))+y[3](y[6](t)));

eq[6] := diff(y[6](t), t) = H2*(y[2](t)+y[3](t))/(y[2](y[6](t))+y[3](y[6](t)))

(14)

H1:=1/2+1/2*tanh(100*(t-35));H2:=1/2+1/2*tanh(100*(t-197));

H1 := 1/2+(1/2)*tanh(100*t-3500)

H2 := 1/2+(1/2)*tanh(100*t-19700)

(15)

alpha:=1.8;beta:=20.;gamma1:=0.002;r:=5.*10^4;s:=10.^5;

alpha := 1.8

beta := 20.

gamma1 := 0.2e-2

r := 50000.

s := 100000.

(16)

seq(eq[i],i=1..6);

diff(y[1](t), t) = -50000.*y[1](t)*y[2](t)-100000.*y[1](t)*y[4](t), diff(y[2](t), t) = -50000.*y[1](t)*y[2](t)+90000.0*y[1](y[5](t))*y[2](y[5](t))*(1/2+(1/2)*tanh(100*t-3500)), diff(y[3](t), t) = 50000.*y[1](t)*y[2](t), diff(y[4](t), t) = -100000.*y[1](t)*y[4](t)-0.2e-2*y[4](t)+1000000.*y[1](y[6](t))*y[2](y[6](t))*(1/2+(1/2)*tanh(100*t-19700)), diff(y[5](t), t) = (1/2+(1/2)*tanh(100*t-3500))*(y[1](t)*y[2](t)+y[3](t))/(y[1](y[5](t))*y[2](y[5](t))+y[3](y[5](t))), diff(y[6](t), t) = (1/2+(1/2)*tanh(100*t-19700))*(y[2](t)+y[3](t))/(y[2](y[6](t))+y[3](y[6](t)))

(17)

ics:=y[1](0)=5.*10^(-6),y[2](0)=10.^(-15),y[3](0)=0,y[4](0)=0,y[5](0)=1e-40,y[6](0)=1e-20;

ics := y[1](0) = 0.5000000000e-5, y[2](0) = 0.1000000000e-14, y[3](0) = 0, y[4](0) = 0, y[5](0) = 0.1e-39, y[6](0) = 0.1e-19

(18)

#infolevel[all]:=10;

sol:=dsolve({seq(eq[i],i=1..6),ics},type=numeric,delaymax=300,initstep=1e-6,abserr=[1e-21,1e-21,1e-21,1e-21,1e-9,1e-9],[y[1](t),y[2](t),y[3](t),y[4](t),y[5](t),y[6](t)],relerr=1e-9,maxstep=10,optimize=false,compile=true,maxfun=0):

 

 

 note that compile = true was used for efficiency

t11:=time():sol(300);time()-t11;

[t = 300., y[1](t) = 0.615611371327094e-15, y[2](t) = 0.337706811581908e-6, y[3](t) = 0.422136411682798e-6, y[4](t) = 0.214253771204037e-5, y[5](t) = 299.999986716780, y[6](t) = 299.643054284209]

.141

(19)

with(plots):

nd:=nops(radar5data);

nd := 1059

(20)

radar5data[nd];

["X", "=", 300.000000, "Y", "=", 0.6154486288e-15, 0.3377120916e-6, 0.4221403310e-6, 0.2142554563e-5, 299.9999999, 299.6430338]

(21)

 Values at t = 300 match with expected results.

pr[1]:=plot([seq([radar5data[i][3],log(radar5data[i][6])/log(10)],i=1..nd)],style=point,color=green):

p[1]:=odeplot(sol,[t,log(y[1](t))/log(10)],0..300,axes=boxed,thickness=3):

display({pr[1],p[1]});

 

pr[2]:=plot([seq([radar5data[i][3],log(radar5data[i][7])/log(10)],i=1..nd)],style=point,color=green):

p[2]:=odeplot(sol,[t,log(y[2](t))/log(10)],0..300,axes=boxed,thickness=3,numpoints=1000):

display({pr[2],p[2]});

 

pr[3]:=plot([seq([radar5data[i][3],log(radar5data[i][8])/log(10)],i=2..nd)],style=point,color=green):

 

p[3]:=odeplot(sol,[t,log(y[3](t))/log(10)],0..300,axes=boxed,thickness=3):

display({pr[3],p[3]});

 

pr[4]:=plot([seq([radar5data[i][3],log(radar5data[i][9])/log(10)],i=496..nd)],style=point,color=green,view=[197..300,-9..-5]):

p[4]:=odeplot(sol,[t,log(y[4](t))/log(10)],197..300,axes=boxed,thickness=3,view=[197..300,-9..-5]):

display({pr[4],p[4]});

 

pr[5]:=plot([seq([radar5data[i][3],radar5data[i][10]],i=1..nd)],style=point,color=green):

p[5]:=odeplot(sol,[t,y[5](t)],0..300,axes=boxed,thickness=3):

display({pr[5],p[5]});

 

pr[6]:=plot([seq([radar5data[i][3],radar5data[i][11]],i=1..nd)],style=point,color=green):

p[6]:=odeplot(sol,[t,y[6](t)],0..300,axes=boxed,thickness=3):

display({pr[6],p[6]});

 


Download delayimmunetopost.mws

Presentations of the first national congress of civil engineering developed at the University Cesar Vallejo. From 10 to 12 November 2014.

 

CONIC_UCV.pdf

(in spanish)

Lenin Araujo Castillo

Physics Pure

Computer Science

 

 

 

 

We have just released a new version of MapleSim.

MapleSim 7 makes it substantially easier to explore and validate designs, create and manage libraries of custom components, and use your MapleSim models with other tools. It includes:

  • Easy model investigation. A new Results Manager gives you greater flexibility when it comes to investigating your simulation results, including the ability to compare simulation runs on the same axes, instantly plot both probed and unprobed variables, and easily create custom plots.
  • Convenient library creation. With MapleSim 7, it is significantly easier to create, manage, and share libraries of custom components.
  • Improved Modelica support. MapleSim 7 expands the support of the Modelica language so that more Modelica definitions can be used directly inside MapleSim.

We have also updated and expanded the MapleSim 7 family of add-on products:

  • The new MapleSim Battery Library, which is available as a separate add-on, allows you to incorporate physics-based predictive models of battery cells into your system models so you can take battery behavior into account early in the design process. 
  • The MapleSim Connector for FMI, which allows engineers to share very efficient, high-fidelity models created in MapleSim with other modeling tools, has been expanded to support more export formats for co-simulation and model exchange.

See What’s New in MapleSim 7 for more information about these and other improvements in MapleSim.

 

eithne

Maplesoft is holding its first ever Virtual User Summit on Feb. 27.  You’ll be able to watch presentations by both Maplesoft and Maplesoft customers, ask questions, have discussions in the lounge with other attendees, and even enter a draw, all from the comfort of your own home or office.

Here’s the agenda.  We’ll release more detailed information on speakers and session times in the next couple of weeks.

For more information and to register:  Maplesoft Virtual User Conference

We're looking forward to seeing you there. (Well, "seeing you" :-))

eithne

Some errors in MapleSim Tutorials! While I was working with MapleSim tutorials, I have encountered some problems that the results did not match the answer. I have uploaded the files, that I have worked on, in a zip file to this address (http://www.mediafire.com/download/2pirvfbawyjfa18/MapleSim_Projects.zip). The problems are as follows:

1. Based on video tutorial on this page (

We have just released a free collection of Maple and MapleSim classroom materials that helps bring modern technology to any introductory course in control design. This collection includes interactive classroom demonstrations that illustrate key concepts; lab projects and assignment questions; and example models ready to be explored.  It’s been designed primarily for instructors, but students should find much of it useful, too.

1 Introduction

Three tanks are connected with two pipes. Each tank is initially filled to a different level. A valve in each pipe opens, and the liquid levels gradually reach equilibrium. Here, we model the system in MapleSim (including the influence of flow inertia), and also derive and solve the analytical equations in Maple.

On Monday, August 6 at 1:31 a.m. EDT, NASA will attempt the landing of a new planetary rover, named Curiosity, on the surface of Mars.  The Mars Science Laboratory project is managed by the NASA Jet Propulsion Laboratory (JPL) in Pasadena, California, a world-renowned center for robotic space exploration and advanced science and engineering.  JPL recently began a widespread adoption of Maplesoft technology, and Maplesoft’s products are expected to help JPL save...

Dr. Gilbert Lai is a mentor for the FIRST Robotics team SWAT 771. He is helping an all girls team from grades 7-12 design a basketball-shooting robot for this year’s annual FIRST Robotics Competition. Dr. Lai is using MapleSim and Maple to help the team understand the principles involved and design their robot. This blog post is part of a series that chronicles the progress of the team.  Posts in the series include:

  • Part 1 - ...

Dr. Gilbert Lai is a mentor for the FIRST Robotics team SWAT 771. He is helping an all girls team from grades 7-12 design a basketball-shooting robot for this year’s annual FIRST Robotics Competition. Dr. Lai is using MapleSim and Maple to help the team understand the principles involved and design their robot. This blog post is part of a series that chronicles the progress of the team.  Posts in the series include:

  • Part 1 - 

Dr. Gilbert Lai is a mentor for the FIRST Robotics team SWAT 771. He is helping an all girls team from grades 7-12 design a basketball-shooting robot for this year’s annual FIRST Robotics Competition. Dr. Lai is using MapleSim and Maple to help the team understand the principles involved and design their robot. This blog post is part of a series that chronicles the progress of the team.  Posts in the series include:

  • Part 1 - 

Liquid flowing in a pipeline has inertia.  If a valve at the end of the pipeline suddenly closes, a pressure surge hits the valve, and travels through the pipeline at the speed of sound. The damping effect of fluid friction gradually attenuates the pressure wave.

This phenomenon is called water hammer and can cause damage significant damage, sometimes even rupturing the pipeline.

The pressure wave often produces audible sound. If you’ve ever heard...

> restart; with(LinearAlgebra); assume(omega, real, omega > 0);
> G := 9;
> z := (xi^2+xi/(1+xi^2))/(1+xi^2);
`output redirected...`> print(); # input placeholder
> C := `<,>`(1-z, seq(sin((n-1)*Pi*z), n = 2 .. G));
`output redirected...`> print(); # input placeholder
> g := Transpose(C);
`output redirected...`> print(); # input placeholder
> A := Multiply(C, g);
`output redirected...`> print(); # input placeholder

> restart; with(LinearAlgebra); assume(omega, real, omega > 0);
> G := 9;
> z := (xi^2+xi/(1+xi^2))/(1+xi^2);
`output redirected...`> print(); # input placeholder
> C := `<,>`(1-z, seq(sin((n-1)*Pi*z), n = 2 .. G));
`output redirected...`> print(); # input placeholder
> g := Transpose(C);
`output redirected...`> print(); # input placeholder
> A := Multiply(C, g);
`output redirected...`> print(); # input placeholder

A prospective customer recently asked if we had a MapleSim model of a double pipe heat exchanger. Heat exchangers are a critical unit operation in the process industries, and accurate models are needed for process control studies.  I couldn't find an appropriate model so I decided to derive the dynamic equations, and implement them using MapleSim's custom component interface.  I'll outline my modeling strategy in this blog post.

1 2 3 Page 2 of 3