Services & Resources / Wolfram Forums / MathGroup Archive
-----

MathGroup Archive 2008

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: NDSolve Problem

  • To: mathgroup at smc.vnet.net
  • Subject: [mg91835] Re: NDSolve Problem
  • From: dh <dh at metrohm.ch>
  • Date: Tue, 9 Sep 2008 06:59:15 -0400 (EDT)
  • References: <g9r4d5$crm$1@smc.vnet.net>


Hi,

we all know that a numerical solution is not 100% correct. Know, we must 

distinguish between absolute and relative error. E.g. An error of 1 may 

be negligable if the function value is 10^6, but hugh if the function is 

1. Therefore, if you plot relative errors, you will see that NDSolve is 

not too bad:

Res={(-x1'[t]+x2[t])/x2[t],(-x2'[t]+c^2 x1[t]-(c^2+p^2) Sin[p t])/x2'[t]};

Plot[Evaluate[Res/.sol/.dp],{t,0.0,2}]



Concerning your second question, where you get something like 10^-16 

instead of zero, this is the old problem of cancelation in sums. As soon 

as you calulate with machine precision your numbers are only 

approximative. If you then subtract 2 large numbers that are only 

approximately the same, you will not get zero.

hope this will help, Daniel



pratip wrote:

> Lets take this system of ODE and solve with DSolve

> 

> exsol = DSolve[{x1'[t] == x2[t], 

>     x2'[t] == c^2 x1[t] - (c^2 + p^2) Sin[p t], x1[0] == 0, 

>     x2[0] == Pi}, {x1[t], x2[t]}, t] // FullSimplify

> 

> solution is as follows

> {{x1[t] -> Sin[p t] + ((-p + \[Pi]) Sinh[c t])/c, 

>   x2[t] -> p Cos[p t] + (-p + \[Pi]) Cosh[c t]}}

> 

> Here is the values of  the parameters

> 

> dp = {c -> 100, p -> 2}

> 

> Now solve with NDSolve

> 

> sol = NDSolve[{x1'[t] == x2[t], 

>     x2'[t] == c^2 x1[t] - (c^2 + p^2) Sin[p t], x1[0] == 0, 

>     x2[0] == Pi} /. dp, {x1, x2}, {t, 0, 2}]

> 

> Lets Plot the residual error

> Res = {-x1'[t] + x2[t], -x2'[t] + c^2 x1[t] - (c^2 + p^2) Sin[p t]};

> Plot[Evaluate[Res /. sol /. dp], {t, 0, 2}]

> 

> One can see the exponential increase in error. Please tell me if it is possible to manage this error for any parameter values. I want to know how I can use NDSolve to simulate such ODE with very little residual error just like any other normal system of ODE. Should I choose some special method options?

> 

> Another very funny thing

> 

> der = D[exsol, t];

> Res /. exsol /. der /. dp // Flatten // FullSimplify

> 

> gives {0,0}. But lets take the residual with out simplifying.

> res = Res /. exsol /. der /. dp // Flatten;

> Plot[Evaluate[res], {t, 0, 1.2}, PlotStyle -> Thick]

> 

> The Plot shows how badly exponential error shows up in the boundary. Is it a problem with Plot command? why it can not understand without simplifying that "res={0,0}". I hope you guys have an answer.

> 





-- 



Daniel Huber

Metrohm Ltd.

Oberdorfstr. 68

CH-9100 Herisau

Tel. +41 71 353 8585, Fax +41 71 353 8907

E-Mail:<mailto:dh at metrohm.com>

Internet:<http://www.metrohm.com>




  • Prev by Date: Re: How can I do a "grep -v" equivalent in Import[]?
  • Next by Date: Re: matrix matching
  • Previous by thread: NDSolve Problem
  • Next by thread: matrix matching